Merge pull request #3868 from hakonhagland/gl_wtest6

Well testing of gas lift wells

opm/simulators/wells/GasLiftSingleWellGeneric.cpp

@@ -294,19 +294,38 @@ checkThpControl_() const
     return thp_control;
 }
 
+std::pair<std::optional<double>,double>
+GasLiftSingleWellGeneric::
+computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const
+{
+    auto alq = this->orig_alq_;
+    double new_alq = alq;
+    std::optional<double> bhp;
+    while (alq <= (this->max_alq_ + this->increment_)) {
+        if (bhp = computeBhpAtThpLimit_(alq); bhp) {
+            new_alq = alq;
+            break;
+        }
+        alq += this->increment_;
+    }
+    return {bhp, new_alq};
+}
 
-std::optional<GasLiftSingleWellGeneric::BasicRates>
+std::pair<std::optional<GasLiftSingleWellGeneric::BasicRates>,double>
 GasLiftSingleWellGeneric::
 computeInitialWellRates_() const
 {
     std::optional<BasicRates> rates;
-    if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) {
+    double initial_alq = this->orig_alq_;
+    //auto alq = initial_alq;
+    //if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) {
+    if (auto [bhp, alq] = computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_(); bhp) {
         {
             const std::string msg = fmt::format(
-                "computed initial bhp {} given thp limit and given alq {}",
+                "computed initial bhp {} given thp limit and given alq {}", *bhp, alq);
-                *bhp, this->orig_alq_);
             displayDebugMessage_(msg);
         }
+        initial_alq = alq;
         auto [new_bhp, bhp_is_limited] = getBhpWithLimit_(*bhp);
         rates = computeWellRates_(new_bhp, bhp_is_limited);
         if (rates) {
@@ -320,7 +339,7 @@ computeInitialWellRates_() const
     else {
         displayDebugMessage_("Aborting optimization.");
     }
-    return rates;
+    return {rates, initial_alq};
 }
 
 std::optional<GasLiftSingleWellGeneric::LimitedRates>
@@ -819,12 +838,13 @@ getRateWithGroupLimit_(
 }
 
 
-std::optional<GasLiftSingleWellGeneric::LimitedRates>
+std::pair<std::optional<GasLiftSingleWellGeneric::LimitedRates>, double>
 GasLiftSingleWellGeneric::
 getInitialRatesWithLimit_() const
 {
     std::optional<LimitedRates> limited_rates;
-    if (auto rates = computeInitialWellRates_(); rates) {
+    double initial_alq = this->orig_alq_;
+    if (auto [rates, alq] = computeInitialWellRates_(); rates) {
         if (this->debug) {
             displayDebugMessage_(
                              "Maybe limiting initial rates before optimize loop..");
@@ -832,8 +852,9 @@ getInitialRatesWithLimit_() const
         auto temp_rates = getLimitedRatesFromRates_(*rates);
         BasicRates old_rates = getWellStateRates_();
         limited_rates = updateRatesToGroupLimits_(old_rates, temp_rates);
+        initial_alq = alq;
     }
-    return limited_rates;
+    return {limited_rates, initial_alq};
 }
 
 GasLiftSingleWellGeneric::LimitedRates
@@ -1166,18 +1187,17 @@ runOptimizeLoop_(bool increase)
 {
     if (this->debug) debugShowProducerControlMode();
     std::unique_ptr<GasLiftWellState> ret_value; // nullptr initially
-    auto rates = getInitialRatesWithLimit_();
+    auto [rates, cur_alq] = getInitialRatesWithLimit_();
     if (!rates) return ret_value;
     // if (this->debug) debugShowBhpAlqTable_();
     if (this->debug) debugShowAlqIncreaseDecreaseCounts_();
     if (this->debug) debugShowTargets_();
     bool success = false;  // did we succeed to increase alq?
     bool alq_is_limited = false;
-    auto cur_alq = this->orig_alq_;
     LimitedRates new_rates = *rates;
     auto [temp_rates2, new_alq] = maybeAdjustALQbeforeOptimizeLoop_(
                                                  *rates, cur_alq, increase);
-    if (checkInitialALQmodified_(new_alq, cur_alq)) {
+    if (checkInitialALQmodified_(new_alq, this->orig_alq_)) {
         auto delta_alq = new_alq - cur_alq;
         new_rates = temp_rates2;
         cur_alq = new_alq;
@@ -1591,7 +1611,10 @@ checkAlqOutsideLimits(double alq, [[maybe_unused]] double oil_rate)
                 // NOTE: checking for an upper limit should not be necessary
                 // when decreasing alq.. so this is just to catch an
                 // illegal state at an early point.
-                if (alq >= this->parent.max_alq_) {
+                if (this->parent.checkALQequal_(alq, this->parent.max_alq_)) {
+                    return false;
+                }
+                else if (alq > this->parent.max_alq_) {
                     warn_( "unexpected: alq above upper limit when trying to "
                         "decrease lift gas. aborting iteration.");
                     result = true;

opm/simulators/wells/GasLiftSingleWellGeneric.hpp

@@ -250,7 +250,8 @@ protected:
     bool checkInitialALQmodified_(double alq, double initial_alq) const;
     bool checkThpControl_() const;
     virtual std::optional<double> computeBhpAtThpLimit_(double alq) const = 0;
-    std::optional<BasicRates> computeInitialWellRates_() const;
+    std::pair<std::optional<double>,double> computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const;
+    std::pair<std::optional<BasicRates>,double> computeInitialWellRates_() const;
     std::optional<LimitedRates> computeLimitedWellRatesWithALQ_(double alq) const;
     virtual BasicRates computeWellRates_(double bhp, bool bhp_is_limited,                                                             bool debug_output = true) const = 0;
     std::optional<BasicRates> computeWellRatesWithALQ_(double alq) const;
@@ -272,7 +273,7 @@ protected:
                            const BasicRates& rates) const;
     std::pair<double, bool> getGasRateWithGroupLimit_(
                            double new_gas_rate, double gas_rate) const;
-    std::optional<LimitedRates> getInitialRatesWithLimit_() const;
+    std::pair<std::optional<LimitedRates>,double> getInitialRatesWithLimit_() const;
     LimitedRates getLimitedRatesFromRates_(const BasicRates& rates) const;
     std::tuple<double,double,bool,bool> getLiquidRateWithGroupLimit_(
                            const double new_oil_rate, const double oil_rate,

opm/simulators/wells/WellInterface.hpp

@@ -248,6 +248,10 @@ public:
 
     void checkWellOperability(const Simulator& ebos_simulator, const WellState& well_state, DeferredLogger& deferred_logger);
 
+    void gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
+                                             WellState& well_state,
+                                             DeferredLogger& deferred_logger);
+
     // check whether the well is operable under the current reservoir condition
     // mostly related to BHP limit and THP limit
     void updateWellOperability(const Simulator& ebos_simulator,

opm/simulators/wells/WellInterface_impl.hpp

@@ -332,6 +332,9 @@ namespace Opm
                 return;
             }
 
+            if (this->isProducer()) {
+                gliftBeginTimeStepWellTestUpdateALQ(simulator, well_state_copy, deferred_logger);
+            }
             updateWellOperability(simulator, well_state_copy, deferred_logger);
             if ( !this->isOperableAndSolvable() ) {
                 const auto msg = fmt::format("WTEST: Well {} is not operable (physical)", this->name());
@@ -585,7 +588,49 @@ namespace Opm
         updateWellOperability(ebos_simulator, well_state, deferred_logger);
     }
 
-
+    template<typename TypeTag>
+    void
+    WellInterface<TypeTag>::
+    gliftBeginTimeStepWellTestUpdateALQ(const Simulator& ebos_simulator,
+                          WellState& well_state,
+                          DeferredLogger& deferred_logger)
+    {
+        const auto& summary_state = ebos_simulator.vanguard().summaryState();
+        const auto& well_name = this->name();
+        if (!this->wellHasTHPConstraints(summary_state)) {
+            const std::string msg = fmt::format("GLIFT WTEST: Well {} does not have THP constraints", well_name);
+            deferred_logger.info(msg);
+            return;
+        }
+        const auto& well_ecl = this->wellEcl();
+        const auto& schedule = ebos_simulator.vanguard().schedule();
+        auto report_step_idx = ebos_simulator.episodeIndex();
+        const auto& glo = schedule.glo(report_step_idx);
+        if (!glo.has_well(well_name)) {
+            const std::string msg = fmt::format(
+                "GLIFT WTEST: Well {} : Gas Lift not activated: "
+                "WLIFTOPT is probably missing. Skipping.", well_name);
+            deferred_logger.info(msg);
+            return;
+        }
+        const auto& gl_well = glo.well(well_name);
+        auto& max_alq_optional = gl_well.max_rate();
+        double max_alq;
+        if (max_alq_optional) {
+            max_alq = *max_alq_optional;
+        }
+        else {
+            const auto& controls = well_ecl.productionControls(summary_state);
+            const auto& table = this->vfpProperties()->getProd()->getTable(controls.vfp_table_number);
+            const auto& alq_values = table.getALQAxis();
+            max_alq = alq_values.back();
+        }
+        well_state.setALQ(well_name, max_alq);
+        const std::string msg = fmt::format(
+            "GLIFT WTEST: Well {} : Setting ALQ to max value: {}",
+            well_name, max_alq);
+        deferred_logger.info(msg);
+    }
 
     template<typename TypeTag>
     void