Merge pull request #6238 from bska/assign-shut-conns-for-owned-wells-only

Assign Shut Connection Data for Owned Wells Only

opm/simulators/flow/GenericOutputBlackoilModule.cpp

@@ -157,7 +157,7 @@ GenericOutputBlackoilModule(const EclipseState& eclState,
     , flowsC_(schedule, summaryConfig)
     , rftC_(eclState_, schedule_,
             [this](const std::string& wname)
-            { return !isDefunctParallelWell(wname); })
+            { return this->isOwnedByCurrentRank(wname); })
     , rst_conv_(std::move(globalCell), comm)
     , local_data_valid_(false)
 {

opm/simulators/flow/GenericOutputBlackoilModule.hpp

@@ -330,7 +330,8 @@ protected:
     void setupExtraBlockData(const std::size_t        reportStepNum,
                              std::function<bool(int)> isCartIdxOnThisRank);
 
-    virtual bool isDefunctParallelWell(std::string wname) const = 0;
+    virtual bool isDefunctParallelWell(const std::string& wname) const = 0;
+    virtual bool isOwnedByCurrentRank(const std::string& wname) const = 0;
 
     const EclipseState& eclState_;
     const Schedule& schedule_;

opm/simulators/flow/OutputBlackoilModule.hpp

@@ -597,7 +597,7 @@ private:
         Problem, std::void_t<decltype(std::declval<Problem>().geoMechModel())>
     > : public std::true_type {};
 
-    bool isDefunctParallelWell(std::string wname) const override
+    bool isDefunctParallelWell(const std::string& wname) const override
     {
         if (simulator_.gridView().comm().size() == 1)
             return false;
@@ -607,6 +607,11 @@ private:
         return candidate == parallelWells.end() || *candidate != value;
     }
 
+    bool isOwnedByCurrentRank(const std::string& wname) const override
+    {
+        return this->simulator_.problem().wellModel().isOwner(wname);
+    }
+
     void updateFluidInPlace_(const ElementContext& elemCtx, const unsigned dofIdx)
     {
         const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, /*timeIdx=*/0);

opm/simulators/flow/OutputCompositionalModule.hpp

@@ -354,7 +354,7 @@ public:
     }
 
 private:
-    bool isDefunctParallelWell(std::string wname) const override
+    bool isDefunctParallelWell(const std::string& wname) const override
     {
         if (simulator_.gridView().comm().size() == 1)
             return false;
@@ -364,6 +364,14 @@ private:
         return candidate == parallelWells.end() || *candidate != value;
     }
 
+    bool isOwnedByCurrentRank(const std::string& wname) const override
+    {
+        // Note: This statement is not correct for distributed wells and
+        // will need additional logic once those are supported for
+        // compositional flows.
+        return ! this->isDefunctParallelWell(wname);
+    }
+
     void createLocalRegion_(std::vector<int>& region)
     {
         std::size_t elemIdx = 0;

opm/simulators/flow/RFTContainer.cpp

@@ -85,19 +85,21 @@ addToWells(data::Wells& wellDatas,
         gatherAndUpdateMap(gasConnectionSaturations_, comm);
     }
 
-    for (const auto& well: schedule_.getWells(reportStepNum)) {
-
-        // don't bother with wells not on this process
-        if (!wellQuery_(well.name())) {
+    for (const auto& wname : this->schedule_.wellNames(reportStepNum)) {
+        // Don't bother with wells not on this process.
+        if (!wellQuery_(wname)) {
             continue;
         }
 
-        //add data infrastructure for shut wells
-        if (!wellDatas.count(well.name())) {
-            auto& wellData = wellDatas[well.name()];
-            wellData.connections.reserve(well.getConnections().size());
-            std::transform(well.getConnections().begin(),
-                           well.getConnections().end(),
+        // Add data infrastructure for shut wells.
+        if (!wellDatas.count(wname)) {
+            const auto& conns = this->schedule_[reportStepNum]
+                .wells(wname).getConnections();
+
+            auto& wellData = wellDatas[wname];
+            wellData.connections.reserve(conns.size());
+
+            std::transform(conns.begin(), conns.end(),
                            std::back_inserter(wellData.connections),
                            [](const auto& connection)
                            {
@@ -107,7 +109,7 @@ addToWells(data::Wells& wellDatas,
                            });
         }
 
-        data::Well& wellData = wellDatas.at(well.name());
+        data::Well& wellData = wellDatas.at(wname);
 
         auto cond_assign = [](double& dest, unsigned idx, const auto& map)
         {

opm/simulators/wells/BlackoilWellModel.hpp

@@ -220,6 +220,8 @@ template<class Scalar> class WellContributions;
 
                 this->assignWellTargets(wsrpt);
 
+                this->assignDynamicWellStatus(wsrpt, this->reportStepIndex());
+
                 // Assigning (a subset of the) property values in shut
                 // connections should be the last step of wellData().
                 this->assignShutConnections(wsrpt, this->reportStepIndex());

opm/simulators/wells/BlackoilWellModelGeneric.cpp

@@ -1025,18 +1025,31 @@ setWsolvent(const Group& group,
     }
 }
 
+template <class Scalar>
+template <typename LoopBody>
+void BlackoilWellModelGeneric<Scalar>::
+loopOwnedWells(LoopBody&& loopBody) const
+{
+    auto wellIndex = 0 * this->wells_ecl_.size();
+
+    for (const auto& pwInfo : this->local_parallel_well_info_) {
+        if (pwInfo.get().isOwner()) {
+            loopBody(wellIndex, this->wells_ecl_[wellIndex]);
+        }
+
+        ++wellIndex;
+    }
+}
+
 template<class Scalar>
 void BlackoilWellModelGeneric<Scalar>::
 assignWellTargets(data::Wells& wsrpt) const
 {
-    auto pwInfo = this->local_parallel_well_info_.begin();
-
-    for (const auto& well : this->wells_ecl_) {
-        if (! pwInfo++->get().isOwner()) {
-            continue;
-        }
-
-        // data::Wells is a std::map<>
+    this->loopOwnedWells([this, &wsrpt]
+                         ([[maybe_unused]] const auto wellIndex,
+                          const Well&                 well)
+    {
+        // data::Wells (i.e., 'wsrpt') is a std::map<>
         auto& limits = wsrpt[well.name()].limits;
 
         if (well.isProducer()) {
@@ -1045,7 +1058,7 @@ assignWellTargets(data::Wells& wsrpt) const
         else {
             this->assignInjectionWellTargets(well, limits);
         }
-    }
+    });
 }
 
 template<class Scalar>
@@ -1093,31 +1106,54 @@ assignInjectionWellTargets(const Well& well, data::WellControlLimits& limits) co
     limits.set(rateItem, ctrl.surface_rate);
 }
 
+template<class Scalar>
+void BlackoilWellModelGeneric<Scalar>::
+assignDynamicWellStatus(data::Wells& wsrpt,
+                        const int reportStepIndex) const
+{
+    this->loopOwnedWells([this, reportStepIndex, &wsrpt]
+                         (const auto wellID, const Well& well)
+    {
+        auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<>
+
+        xwel.dynamicStatus = this->schedule()[reportStepIndex]
+            .wells(well.name()).getStatus();
+
+        if ((xwel.dynamicStatus == Well::Status::OPEN) &&
+            this->wellTestState().well_is_closed(well.name()) &&
+            !this->wasDynamicallyShutThisTimeStep(wellID))
+        {
+            // Well is supposed to be flowing according to the run
+            // specification (Schedule object), but it is not operable or
+            // cannot meet its economic limits (well testing).
+            //
+            // Assign status based on the well's defined automatic shut-in
+            // procedure.
+            xwel.dynamicStatus = well.getAutomaticShutIn()
+                ? Well::Status::SHUT : Well::Status::STOP;
+        }
+    });
+}
+
 template<class Scalar>
 void BlackoilWellModelGeneric<Scalar>::
 assignShutConnections(data::Wells& wsrpt,
                       const int reportStepIndex) const
 {
-    auto wellID = 0;
+    this->loopOwnedWells([this, reportStepIndex, &wsrpt]
+                         ([[maybe_unused]] const auto wellID,
+                          const Well&                 well)
+    {
+        const auto wellIsOpen =
+            this->schedule()[reportStepIndex].wells(well.name())
+            .getStatus() == Well::Status::OPEN;
 
-    for (const auto& well : this->wells_ecl_) {
-        auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<>
-
-        xwel.dynamicStatus = this->schedule()
-            .getWell(well.name(), reportStepIndex).getStatus();
-
-        const auto wellIsOpen = xwel.dynamicStatus == Well::Status::OPEN;
         auto skip = [wellIsOpen](const Connection& conn)
         {
             return wellIsOpen && (conn.state() != Connection::State::SHUT);
         };
 
-        if (this->wellTestState().well_is_closed(well.name()) &&
-            !this->wasDynamicallyShutThisTimeStep(wellID))
-        {
-            xwel.dynamicStatus = well.getAutomaticShutIn()
-                ? Well::Status::SHUT : Well::Status::STOP;
-        }
+        auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<>
 
         auto& xcon = xwel.connections;
         for (const auto& conn : well.getConnections()) {
@@ -1133,9 +1169,7 @@ assignShutConnections(data::Wells& wsrpt,
             xc.trans_factor = conn.CF();
             xc.d_factor = conn.dFactor();
         }
-
-        ++wellID;
-    }
+    });
 }
 
 template<class Scalar>

opm/simulators/wells/BlackoilWellModelGeneric.hpp

@@ -36,12 +36,14 @@
 #include <opm/simulators/wells/ConnectionIndexMap.hpp>
 #include <opm/simulators/wells/ParallelPAvgDynamicSourceData.hpp>
 #include <opm/simulators/wells/ParallelWBPCalculation.hpp>
+#include <opm/simulators/wells/ParallelWellInfo.hpp>
 #include <opm/simulators/wells/PerforationData.hpp>
 #include <opm/simulators/wells/WellFilterCake.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/WellTracerRate.hpp>
 #include <opm/simulators/wells/WGState.hpp>
 
+#include <algorithm>
 #include <cstddef>
 #include <functional>
 #include <map>
@@ -61,7 +63,6 @@ namespace Opm {
     template<class Scalar> class GasLiftWellState;
     class Group;
     class GuideRateConfig;
-    template<class Scalar> class ParallelWellInfo;
     class RestartValue;
     class Schedule;
     struct SimulatorUpdate;
@@ -268,6 +269,17 @@ public:
     parallelWellInfo(const std::size_t idx) const
     { return local_parallel_well_info_[idx].get(); }
 
+    bool isOwner(const std::string& wname) const
+    {
+        auto pwInfoPos = std::find_if(this->parallel_well_info_.begin(),
+                                      this->parallel_well_info_.end(),
+                                      [&wname](const auto& pwInfo)
+                                      { return pwInfo.name() == wname; });
+
+        return (pwInfoPos != this->parallel_well_info_.end())
+            && pwInfoPos->isOwner();
+    }
+
     const ConnectionIndexMap& connectionIndexMap(const std::size_t idx)
     { return conn_idx_map_[idx]; }
 
@@ -370,11 +382,33 @@ protected:
                                   const int pvtreg,
                                   std::vector<Scalar>& resv_coeff) = 0;
 
+    /// Assign dynamic well status for each well owned by current rank
+    ///
+    /// \param[in,out] wsrpt Well solution object.  On exit, holds current
+    /// values for \code data::Well::dynamicStatus \endcode.
+    ///
+    /// \param[in] reportStepIdx Zero-based index of current report step.
+    void assignDynamicWellStatus(data::Wells& wsrpt,
+                                 const int reportStepIdx) const;
+
+    /// Assign basic result quantities for shut connections of wells owned
+    /// by current rank.
+    ///
+    /// Mostly provided for summary file output purposes.  Applies to fully
+    /// shut/stopped wells and shut connections of open/flowing wells.
+    ///
+    /// \param[in,out] wsrpt Well solution object.  On exit, also contains a
+    /// few quantities, like the D factor, the Kh product and the CTF, for
+    /// shut connections.
+    ///
+    /// \param[in] reportStepIdx Zero-based index of current report step.
     void assignShutConnections(data::Wells& wsrpt,
                                const int reportStepIndex) const;
+
     void assignWellTargets(data::Wells& wsrpt) const;
     void assignProductionWellTargets(const Well& well, data::WellControlLimits& limits) const;
     void assignInjectionWellTargets(const Well& well, data::WellControlLimits& limits) const;
+
     void assignGroupControl(const Group& group,
                             data::GroupData& gdata) const;
     void assignGroupValues(const int reportStepIdx,
@@ -550,8 +584,25 @@ private:
 
     void updateEclWellsCTFFromAction(const int timeStepIdx,
                                      const SimulatorUpdate& sim_update);
-};
 
+    /// Run caller-defined code for each well owned by current rank
+    ///
+    /// 'const' version.
+    ///
+    /// \tparam LoopBody Call-back type for user-defined code.  Expected to
+    /// support a function call operator of the form
+    /// \code
+    ///   void operator()(const std::size_t i, const Well& well) const
+    /// \endcode
+    /// which will be invoked for each well owned by the local process.  The
+    /// parameters are the index into \c wells_ecl_ and the object at that
+    /// index, respectively.
+    ///
+    /// \param[in] loopBody Call-back function.  Typically defined as a
+    /// lambda expression.
+    template <typename LoopBody>
+    void loopOwnedWells(LoopBody&& loopBody) const;
+};
 
 } // namespace Opm