Record Dynamic Well Status in 'wellData()' Output

This commit sets the 'data::Well::dynamicStatus' based on the
dynamically updated 'Schedule' object (i.e., from ACTIONX and
similar) and the results of well/operability testing (WECON and/or
WTEST).  If a well is closed due to economic limits (WECON) we still
provide summary-style data at the timestep that closed the well, but
omit this data at later steps until the well reopens.

We add a new parameter to WellState::report() to distinguish these
situations.

This is in preparation of making the 'BlackoilWellModel' manage both
open and shut wells alike.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -217,7 +217,12 @@ namespace Opm {
 
             Opm::data::Wells wellData() const
             {
-                auto wsrpt = well_state_.report(phase_usage_, Opm::UgGridHelpers::globalCell(grid()));
+                auto wsrpt = well_state_
+                    .report(phase_usage_, Opm::UgGridHelpers::globalCell(grid()),
+                            [this](const int well_ndex) -> bool
+                    {
+                        return this->wasDynamicallyShutThisTimeStep(well_ndex);
+                    });
 
                 this->assignWellGuideRates(wsrpt);
                 this->assignShutConnections(wsrpt);
@@ -347,6 +352,7 @@ namespace Opm {
             std::unique_ptr<GuideRate> guideRate_;
 
             std::map<std::string, double> node_pressures_{}; // Storing network pressures for output.
+            mutable std::unordered_set<std::string> closed_this_step_{};
 
             // used to better efficiency of calcuation
             mutable BVector scaleAddRes_;
@@ -430,6 +436,8 @@ namespace Opm {
 
             int numPhases() const;
 
+            int reportStepIndex() const;
+
             void assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger);
 
             // some preparation work, mostly related to group control and RESV,
@@ -489,6 +497,8 @@ namespace Opm {
 
             void runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger);
 
+            bool wasDynamicallyShutThisTimeStep(const int well_index) const;
+
             void assignWellGuideRates(data::Wells& wsrpt) const;
             void assignShutConnections(data::Wells& wsrpt) const;
             void assignGroupValues(const int                               reportStepIdx,

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -504,6 +504,7 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     timeStepSucceeded(const double& simulationTime, const double dt)
     {
+        this->closed_this_step_.clear();
 
         // time step is finished and we are not any more at the beginning of an report step
         report_step_starts_ = false;
@@ -761,6 +762,20 @@ namespace Opm {
             for (int w = 0; w < nw; ++w) {
                 const Well& well_ecl = wells_ecl_[w];
                 const std::string& well_name = well_ecl.name();
+                const auto well_status = this->schedule()
+                    .getWell(well_name, time_step).getStatus();
+
+                if ((well_ecl.getStatus() == Well::Status::SHUT) ||
+                    (well_status          == Well::Status::SHUT))
+                {
+                    // Due to ACTIONX the well might have been closed behind our back.
+                    if (well_ecl.getStatus() != Well::Status::SHUT) {
+                        this->closed_this_step_.insert(well_name);
+                        well_state_.shutWell(w);
+                    }
+
+                    continue;
+                }
 
                 // A new WCON keywords can re-open a well that was closed/shut due to Physical limit
                 if (this->wellTestState_.hasWellClosed(well_name)) {
@@ -796,13 +811,6 @@ namespace Opm {
                     }
                 }
 
-                // Due to ACTIONX the well might have been closed 'behind our back'.
-                const auto well_status = schedule().getWell(well_name, time_step).getStatus();
-                if (well_status == Well::Status::SHUT) {
-                    well_state_.shutWell(w);
-                    continue;
-                }
-
                 // If a production well disallows crossflow and its
                 // (prediction type) rate control is zero, then it is effectively shut.
                 if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) {
@@ -1410,8 +1418,15 @@ namespace Opm {
     {
         Opm::DeferredLogger local_deferredLogger;
         for (const auto& well : well_container_) {
+            const auto wasClosed = wellTestState.hasWellClosed(well->name());
+
             well->updateWellTestState(well_state_, simulationTime, /*writeMessageToOPMLog=*/ true, wellTestState, local_deferredLogger);
+
+            if (!wasClosed && wellTestState.hasWellClosed(well->name())) {
+                this->closed_this_step_.insert(well->name());
+            }
         }
+
         Opm::DeferredLogger global_deferredLogger = gatherDeferredLogger(local_deferredLogger);
         if (terminal_output_) {
             global_deferredLogger.logMessages();
@@ -2680,6 +2695,21 @@ namespace Opm {
     }
 
 
+
+
+
+    template <typename TypeTag>
+    int
+    BlackoilWellModel<TypeTag>::
+    reportStepIndex() const
+    {
+        return std::max(this->ebosSimulator_.episodeIndex(), 0);
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -2746,6 +2776,19 @@ namespace Opm {
 
 
 
+    template <typename TypeTag>
+    bool
+    BlackoilWellModel<TypeTag>::
+    wasDynamicallyShutThisTimeStep(const int well_index) const
+    {
+        return this->closed_this_step_.find(this->wells_ecl_[well_index].name())
+            != this->closed_this_step_.end();
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -2826,15 +2869,30 @@ namespace Opm {
     BlackoilWellModel<TypeTag>::
     assignShutConnections(data::Wells& wsrpt) const
     {
+        auto wellID = 0;
+
         for (const auto& well : this->wells_ecl_) {
-            auto xwPos = wsrpt.find(well.name());
-            if (xwPos == wsrpt.end()) { // No well results.  Unexpected.
-                continue;
+            auto& xwel = wsrpt[well.name()]; // data::Wells is a std::map<>
+
+            xwel.dynamicStatus = this->schedule()
+                .getWell(well.name(), this->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_.hasWellClosed(well.name()) &&
+                !this->wasDynamicallyShutThisTimeStep(wellID))
+            {
+                xwel.dynamicStatus = well.getAutomaticShutIn()
+                    ? Well::Status::SHUT : Well::Status::STOP;
             }
 
-            auto& xcon = xwPos->second.connections;
+            auto& xcon = xwel.connections;
             for (const auto& conn : well.getConnections()) {
-                if (conn.state() != Connection::State::SHUT) {
+                if (skip(conn)) {
                     continue;
                 }
 
@@ -2845,6 +2903,8 @@ namespace Opm {
                 xc.effective_Kh = conn.Kh();
                 xc.trans_factor = conn.CF();
             }
+
+            ++wellID;
         }
     }
 

opm/simulators/wells/WellState.hpp

@@ -30,6 +30,7 @@
 #include <array>
 #include <cassert>
 #include <cstddef>
+#include <functional>
 #include <map>
 #include <memory>
 #include <stdexcept>
@@ -207,6 +208,7 @@ namespace Opm
 
             return well_info.isOwner();
         }
+
         /// The number of wells present.
         int numWells() const
         {
@@ -237,15 +239,21 @@ namespace Opm
             this->thp_[well_index] = 0;
         }
 
-        virtual data::Wells report(const PhaseUsage& pu, const int* globalCellIdxMap) const
+        virtual data::Wells
+        report(const PhaseUsage& pu,
+               const int* globalCellIdxMap,
+               const std::function<bool(const int)>& wasDynamicallyClosed) const
         {
             using rt = data::Rates::opt;
 
             data::Wells dw;
             for( const auto& itr : this->wellMap_ ) {
                 const auto well_index = itr.second[ 0 ];
-                if (this->status_[well_index] == Well::Status::SHUT)
+                if ((this->status_[well_index] == Well::Status::SHUT) &&
+                    ! wasDynamicallyClosed(well_index))
+                {
                     continue;
+                }
 
                 const auto& pwinfo = *parallel_well_info_[well_index];
                 using WellT = std::remove_reference_t<decltype(dw[ itr.first ])>;

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -32,6 +32,7 @@
 #include <algorithm>
 #include <array>
 #include <cassert>
+#include <functional>
 #include <iostream>
 #include <map>
 #include <numeric>
@@ -546,11 +547,13 @@ namespace Opm
             return it->second;
         }
 
-
-
-        data::Wells report(const PhaseUsage &pu, const int* globalCellIdxMap) const override
+        data::Wells
+        report(const PhaseUsage &pu,
+               const int* globalCellIdxMap,
+               const std::function<bool(const int)>& wasDynamicallyClosed) const override
         {
-            data::Wells res = WellState::report(pu, globalCellIdxMap);
+            data::Wells res =
+                WellState::report(pu, globalCellIdxMap, wasDynamicallyClosed);
 
             const int nw = this->numWells();
             if (nw == 0) {
@@ -581,8 +584,9 @@ namespace Opm
 
             for (const auto& wt : this->wellMap()) {
                 const auto w = wt.second[ 0 ];
-                const auto& pwinfo = *parallel_well_info_[w];
-                if ((this->status_[w] == Well::Status::SHUT) || !pwinfo.isOwner())
+                if (((this->status_[w] == Well::Status::SHUT) &&
+                     ! wasDynamicallyClosed(w)) ||
+                    ! this->parallel_well_info_[w]->isOwner())
                 {
                     continue;
                 }

tests/test_wellstatefullyimplicitblackoil.cpp

@@ -274,7 +274,7 @@ BOOST_AUTO_TEST_CASE(Pressure)
 
     setSegPress(wells, wstate);
 
-    const auto rpt = wstate.report(setup.pu, setup.grid.c_grid()->global_cell);
+    const auto rpt = wstate.report(setup.pu, setup.grid.c_grid()->global_cell, [](const int){return false;});
 
     {
         const auto& xw = rpt.at("INJE01");
@@ -323,7 +323,7 @@ BOOST_AUTO_TEST_CASE(Rates)
 
     setSegRates(wells, pu, wstate);
 
-    const auto rpt = wstate.report(pu, setup.grid.c_grid()->global_cell);
+    const auto rpt = wstate.report(pu, setup.grid.c_grid()->global_cell, [](const int){return false;});
 
     const auto wat = pu.phase_used[Opm::BlackoilPhases::Aqua];
     const auto oil = pu.phase_used[Opm::BlackoilPhases::Liquid];