Run WELPI Scaling At Beginning of Report Step

This commit implements the WELPI feature.  We calculate new PI/II
values for all wells in the event of a WELPI request and use those
values for well-specific WELPI request, to calculate CTF scaling
factors.  We then apply those factors to all subsequent editions of
the well provided the connection factors are eligible for
WELPI-based rescaling.

If we trigger a rescaling event we also reset the WellState's
internal copies of the CTFs and reinitialize the Well PI calculators
to ensure the rescaling takes effect immediately.  Since we rely on
PI values being available at the end of each time step we must also
take care to forward those values from the WellState of one report
step to the WellState of the next report step.

Finally, take care not to redo a WELPI scaling if we've already
performed the scaling operation and restart a report step.  This,
in turn, happens if WELPI is requested on the first report step.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -34,6 +34,7 @@
 #include <functional>
 #include <map>
 #include <memory>
+#include <optional>
 #include <string>
 #include <tuple>
 #include <unordered_map>
@@ -322,6 +323,9 @@ namespace Opm {
             bool report_step_starts_;
             bool glift_debug = false;
             bool alternative_well_rate_init_;
+
+            std::optional<int> last_run_wellpi_{};
+
             std::unique_ptr<RateConverterType> rateConverter_;
             std::unique_ptr<VFPProperties<VFPInjProperties,VFPProdProperties>> vfp_properties_;
 
@@ -468,6 +472,8 @@ namespace Opm {
 
             void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent);
 
+            void runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger);
+
             void assignWellGuideRates(data::Wells& wsrpt) const;
             void assignShutConnections(data::Wells& wsrpt) const;
             void assignGroupValues(const int                               reportStepIdx,

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -23,6 +23,8 @@
 #include <opm/simulators/utils/DeferredLoggingErrorHelpers.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 
+#include <opm/parser/eclipse/Units/UnitSystem.hpp>
+
 #include <algorithm>
 #include <utility>
 #include <fmt/format.h>
@@ -238,7 +240,6 @@ namespace Opm {
         // Make wells_ecl_ contain only this partition's non-shut wells.
         wells_ecl_ = getLocalNonshutWells(timeStepIdx, globalNumWells);
 
-        this->initializeWellProdIndCalculators();
         // The well state initialize bhp with the cell pressure in the top cell.
         // We must therefore provide it with updated cell pressures
         this->initializeWellPerfData();
@@ -287,9 +288,13 @@ namespace Opm {
                                    schedule().getVFPInjTables(timeStepIdx),
                                    schedule().getVFPProdTables(timeStepIdx)) );
 
+        this->initializeWellProdIndCalculators();
+        if (this->schedule().getEvents().hasEvent(ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX, timeStepIdx)) {
+            this->runWellPIScaling(timeStepIdx, local_deferredLogger);
+        }
+
         // update the previous well state. This is used to restart failed steps.
         previous_well_state_ = well_state_;
-
     }
 
 
@@ -2550,6 +2555,127 @@ namespace Opm {
 
 
 
+
+
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    runWellPIScaling(const int timeStepIdx, DeferredLogger& local_deferredLogger)
+    {
+        if (this->last_run_wellpi_.has_value() && (*this->last_run_wellpi_ == timeStepIdx)) {
+            // We've already run WELPI scaling for this report step.  Most
+            // common for the very first report step.  Don't redo WELPI scaling.
+            return;
+        }
+
+        auto hasWellPIEvent = [this, timeStepIdx](const int well_index) -> bool
+        {
+            return this->schedule()
+                .hasWellGroupEvent(this->wells_ecl_[well_index].name(),
+                                   ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX,
+                                   timeStepIdx);
+        };
+
+        auto getWellPI = [this](const int well_index) -> double
+        {
+            const auto& pu = this->phase_usage_;
+            const auto  np = this->numPhases();
+            const auto* pi = &this->well_state_.productivityIndex()[np*well_index + 0];
+
+            const auto preferred = this->wells_ecl_[well_index].getPreferredPhase();
+            switch (preferred) { // Should really have LIQUID = OIL + WATER here too...
+            case Phase::WATER:
+                return pu.phase_used[BlackoilPhases::PhaseIndex::Aqua]
+                    ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Aqua]]
+                    : 0.0;
+
+            case Phase::OIL:
+                return pu.phase_used[BlackoilPhases::PhaseIndex::Liquid]
+                    ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Liquid]]
+                    : 0.0;
+
+            case Phase::GAS:
+                return pu.phase_used[BlackoilPhases::PhaseIndex::Vapour]
+                    ? pi[pu.phase_pos[BlackoilPhases::PhaseIndex::Vapour]]
+                    : 0.0;
+
+            default:
+                throw std::invalid_argument {
+                    "Unsupported preferred phase " +
+                    std::to_string(static_cast<int>(preferred))
+                };
+            }
+        };
+
+        auto getWellPIScalingFactor = [this](const int    well_index,
+                                             const double newWellPI) -> double
+        {
+            return this->wells_ecl_[well_index].getWellPIScalingFactor(newWellPI);
+        };
+
+        auto rescaleWellPI =
+            [this, timeStepIdx](const int    well_index,
+                                const double scalingFactor) -> void
+        {
+            {
+                const auto& wname = this->wells_ecl_[well_index].name();
+                auto& schedule = this->ebosSimulator_.vanguard().schedule(); // Mutable
+
+                schedule.applyWellProdIndexScaling(wname, timeStepIdx, scalingFactor);
+                this->wells_ecl_[well_index] = schedule.getWell(wname, timeStepIdx);
+            }
+
+            const auto& well = this->wells_ecl_[well_index];
+            auto& pd     = this->well_perf_data_[well_index];
+            auto  pdIter = pd.begin();
+            for (const auto& conn : well.getConnections()) {
+                if (conn.state() != Connection::State::SHUT) {
+                    pdIter->connection_transmissibility_factor = conn.CF();
+                    ++pdIter;
+                }
+            }
+
+            this->well_state_.resetConnectionTransFactors(well_index, pd);
+            this->prod_index_calc_[well_index].reInit(well);
+        };
+
+        // Minimal well setup to compute PI/II values
+        {
+            auto saved_previous_well_state = this->previous_well_state_;
+            this->previous_well_state_ = this->well_state_;
+
+            well_container_ = createWellContainer(timeStepIdx);
+
+            for (auto& well : well_container_) {
+                well->init(&phase_usage_, depth_, gravity_, local_num_cells_);
+            }
+
+            std::fill(is_cell_perforated_.begin(), is_cell_perforated_.end(), false);
+            for (auto& well : well_container_) {
+                well->updatePerforatedCell(is_cell_perforated_);
+            }
+
+            this->calculateProductivityIndexValues(local_deferredLogger);
+
+            this->previous_well_state_ = std::move(saved_previous_well_state);
+        }
+
+        const auto nw = this->numLocalWells();
+        for (auto wellID = 0*nw; wellID < nw; ++wellID) {
+            if (hasWellPIEvent(wellID)) {
+                const auto newWellPI = getWellPI(wellID);
+                const auto scalingFactor = getWellPIScalingFactor(wellID, newWellPI);
+                rescaleWellPI(wellID, scalingFactor);
+            }
+        }
+
+        this->last_run_wellpi_ = timeStepIdx;
+    }
+
+
+
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -270,6 +270,16 @@ namespace Opm
                                 }
                             }
                         }
+
+                        // Productivity index.
+                        {
+                            auto*       thisWellPI = &this     ->productivityIndex()[newIndex*np + 0];
+                            const auto* thatWellPI = &prevState->productivityIndex()[oldIndex*np + 0];
+
+                            for (int p = 0; p < np; ++p) {
+                                thisWellPI[p] = thatWellPI[p];
+                            }
+                        }
                     }
 
                     // If in the new step, there is no THP related target/limit anymore, its thp value should be