Well testing of gas lift wells

Assign a maximum ALQ value to each GLIFT producer when doing well testing
in beginTimeStep(). This allows the well to be considered open. Then,
later in the timestep, when assemble() is called, the full gas lift
optimization procedure can adjust the ALQ to its correct value.

It is also observed that in some cases when gas lift is switched off by
setting ALQ to zero, and later in the schedule is switched back on again,
it might not be possible to determine bhp from thp for low small ALQ values.
Instead of aborting the gas lift optimization, we should try increasing
ALQ until we get convergence or until the maximum ALQ for the well is
reached.

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