Merge pull request #2824 from hakonhagland/gas_lift3

Implements support for gas lift optimization.

opm/simulators/wells/BlackoilWellModel.hpp

@@ -306,7 +306,7 @@ namespace Opm {
             std::vector<double> depth_;
             bool initial_step_;
             bool report_step_starts_;
-
+            bool glift_debug = false;
             std::unique_ptr<RateConverterType> rateConverter_;
             std::unique_ptr<VFPProperties<VFPInjProperties,VFPProdProperties>> vfp_properties_;
 
@@ -327,6 +327,10 @@ namespace Opm {
             const Schedule& schedule() const
             { return ebosSimulator_.vanguard().schedule(); }
 
+            void gliftDebug(
+                const std::string &msg,
+                Opm::DeferredLogger& deferred_logger) const;
+
             // compute the well fluxes and assemble them in to the reservoir equations as source terms
             // and in the well equations.
             void assemble(const int iterationIdx,

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -26,6 +26,7 @@
 
 #include <utility>
 #include <algorithm>
+#include <fmt/format.h>
 
 namespace Opm {
     template<typename TypeTag>
@@ -330,7 +331,7 @@ namespace Opm {
         Opm::DeferredLogger local_deferredLogger;
 
         well_state_ = previous_well_state_;
-
+        well_state_.disableGliftOptimization();
         const int reportStepIdx = ebosSimulator_.episodeIndex();
         const double simulationTime = ebosSimulator_.time();
 
@@ -393,8 +394,22 @@ namespace Opm {
         //compute well guideRates
         const auto& comm = ebosSimulator_.vanguard().grid().comm();
         WellGroupHelpers::updateGuideRatesForWells(schedule(), phase_usage_, reportStepIdx, simulationTime, well_state_, comm, guideRate_.get());
+        logAndCheckForExceptionsAndThrow(local_deferredLogger,
+            exception_thrown, "beginTimeStep() failed.", terminal_output_);
+
     }
 
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::gliftDebug(
+        const std::string &msg, Opm::DeferredLogger &deferred_logger) const
+    {
+        if (this->glift_debug) {
+            const std::string message = fmt::format(
+                "  GLIFT (DEBUG) : BlackoilWellModel : {}", msg);
+            deferred_logger.info(message);
+        }
+    }
 
     template<typename TypeTag>
     void
@@ -813,6 +828,12 @@ namespace Opm {
              const double dt)
     {
 
+        Opm::DeferredLogger local_deferredLogger;
+        if (this->glift_debug) {
+            const std::string msg = fmt::format(
+                "assemble() : iteration {}" , iterationIdx);
+            gliftDebug(msg, local_deferredLogger);
+        }
         last_report_ = SimulatorReportSingle();
         Dune::Timer perfTimer;
         perfTimer.start();
@@ -821,7 +842,6 @@ namespace Opm {
             return;
         }
 
-        Opm::DeferredLogger local_deferredLogger;
 
         updatePerforationIntensiveQuantities();
 
@@ -855,8 +875,10 @@ namespace Opm {
                 // basically, this is a more updated state from the solveWellEq based on fixed
                 // reservoir state, will tihs be a better place to inialize the explict information?
             }
-
+            gliftDebug("assemble() : running assembleWellEq()..", local_deferredLogger);
+            well_state_.enableGliftOptimization();
             assembleWellEq(B_avg, dt, local_deferredLogger);
+            well_state_.disableGliftOptimization();
 
         } catch (std::exception& e) {
             exception_thrown = 1;
@@ -873,6 +895,8 @@ namespace Opm {
     assembleWellEq(const std::vector<Scalar>& B_avg, const double dt, Opm::DeferredLogger& deferred_logger)
     {
         for (auto& well : well_container_) {
+            well->maybeDoGasLiftOptimization(
+                 well_state_, ebosSimulator_, deferred_logger);
             well->assembleWellEq(ebosSimulator_, B_avg, dt, well_state_, deferred_logger);
         }
     }

opm/simulators/wells/GasLiftRuntime.hpp

@@ -0,0 +1,154 @@
+/*
+  Copyright 2020 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_GASLIFT_RUNTIME_HEADER_INCLUDED
+#define OPM_GASLIFT_RUNTIME_HEADER_INCLUDED
+
+#include <opm/models/utils/propertysystem.hh>
+#include <opm/models/utils/parametersystem.hh>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/output/data/Wells.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+#include <opm/simulators/utils/DeferredLogger.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/GasLiftOpt.hpp>
+// NOTE: StandardWell.hpp includes ourself (GasLiftRuntime.hpp), so we need
+//   to forward declare StandardWell for it to be defined in this file.
+namespace Opm {
+    template<typename TypeTag> class StandardWell;
+}
+#include <opm/simulators/wells/StandardWell.hpp>
+
+#include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+
+#include <cassert>
+#include <iostream>
+#include <map>
+#include <memory>
+#include <optional>
+#include <string>
+#include <vector>
+#include <fmt/format.h>
+
+namespace Opm
+{
+    template<class TypeTag>
+    class GasLiftRuntime {
+        typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
+        using WellState = WellStateFullyImplicitBlackoil;
+        using StdWell = Opm::StandardWell<TypeTag>;
+        // TODO: same definition with WellInterface, and
+        //  WellStateFullyImplicitBlackoil eventually they should go
+        //  to a common header file.
+        static const int Water = BlackoilPhases::Aqua;
+        static const int Oil = BlackoilPhases::Liquid;
+        static const int Gas = BlackoilPhases::Vapour;
+        struct OptimizeState;
+    public:
+        GasLiftRuntime(
+            const StdWell &std_well,
+            const Simulator &ebos_simulator,
+            const SummaryState &summary_state,
+            DeferredLogger &deferred_logger,
+            std::vector<double> &potentials,
+            const WellState &well_state,
+            const Well::ProductionControls &controls
+        );
+        void runOptimize();
+    private:
+        void computeInitialWellRates_();
+        void computeWellRates_(double bhp, std::vector<double> &potentials);
+        void debugShowIterationInfo_(OptimizeState &state, double alq);
+        void debugShowStartIteration_(double alq, bool increase);
+        void displayDebugMessage_(const std::string &msg);
+        void displayWarning_();
+        void displayWarning_(std::string warning);
+        double getGasRateWithLimit_(std::vector<double> &potentials);
+        double getOilRateWithLimit_(std::vector<double> &potentials);
+        void logSuccess_();
+        bool runOptimizeLoop_(bool increase);
+        void setAlqMaxRate_(const GasLiftOpt::Well &well);
+        void setAlqMinRate_(const GasLiftOpt::Well &well);
+        bool tryIncreaseLiftGas_();
+        bool tryDecreaseLiftGas_();
+        void updateWellStateAlqFixedValue_(const GasLiftOpt::Well &well);
+        bool useFixedAlq_(const GasLiftOpt::Well &well);
+        void warnMaxIterationsExceeded_();
+
+        const Well::ProductionControls &controls_;
+        DeferredLogger &deferred_logger_;
+        const Simulator &ebos_simulator_;
+        std::vector<double> &potentials_;
+        const StdWell &std_well_;
+        const SummaryState &summary_state_;
+        const WellState &well_state_;
+        std::string well_name_;
+        bool debug;  // extra debug output
+
+        double alpha_w_;
+        double alpha_g_;
+        double eco_grad_;
+        int gas_pos_;
+        bool has_run_init_ = false;
+        double increment_;
+        double max_alq_;
+        int max_iterations_;
+        double min_alq_;
+        double new_alq_;
+        int oil_pos_;
+        bool optimize_;
+        double orig_alq_;
+        int water_pos_;
+
+        struct OptimizeState {
+            OptimizeState( GasLiftRuntime &parent_, bool increase_ ) :
+                parent(parent_),
+                increase(increase_),
+                it(0),
+                stop_iteration(false),
+                bhp(-1)
+            {}
+
+            GasLiftRuntime &parent;
+            bool increase;
+            int it;
+            bool stop_iteration;
+            double bhp;
+
+            double addOrSubtractAlqIncrement(double alq);
+            double calcGradient(double oil_rate, double new_oil_rate,
+                double gas_rate, double new_gas_rate);
+            bool checkAlqOutsideLimits(double alq, double oil_rate);
+            bool checkEcoGradient(double gradient);
+            bool checkOilRateExceedsTarget(double oil_rate);
+            bool checkRate(double rate, double limit, const std::string rate_str);
+            bool checkWellRatesViolated(std::vector<double> &potentials);
+            bool computeBhpAtThpLimit(double alq);
+            double getBhpWithLimit();
+            void warn_(std::string msg) {parent.displayWarning_(msg);}
+        };
+
+    };
+
+} // namespace Opm
+
+#include "GasLiftRuntime_impl.hpp"
+
+#endif // OPM_GASLIFT_RUNTIME_HEADER_INCLUDED

opm/simulators/wells/GasLiftRuntime_impl.hpp

@@ -0,0 +1,792 @@
+/*
+  Copyright 2020 Equinor ASA.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#include <opm/simulators/wells/StandardWell.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
+#include <opm/simulators/utils/DeferredLogger.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/SummaryState.hpp>
+#include <opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp>
+#include <opm/parser/eclipse/EclipseState/Schedule/GasLiftOpt.hpp>
+
+#include <optional>
+#include <string>
+
+template<typename TypeTag>
+Opm::GasLiftRuntime<TypeTag>::
+GasLiftRuntime(
+    const StdWell &std_well,
+    const Simulator &ebos_simulator,
+    const SummaryState &summary_state,
+    DeferredLogger &deferred_logger,
+    std::vector<double> &potentials,
+    const WellState &well_state,
+    const Well::ProductionControls &controls
+) :
+    controls_{controls},
+    deferred_logger_{deferred_logger},
+    ebos_simulator_{ebos_simulator},
+    potentials_{potentials},
+    std_well_{std_well},
+    summary_state_{summary_state},
+    well_state_{well_state},
+    debug{false}  // extra debugging output
+{
+    int well_index = this->std_well_.indexOfWell();
+    const Well::ProducerCMode& control_mode
+        = well_state_.currentProductionControls()[well_index];
+    if (control_mode != Well::ProducerCMode::THP)
+        throw std::logic_error("Bug in flow - invalid control mode detected\n");
+    const Opm::Schedule& schedule = this->ebos_simulator_.vanguard().schedule();
+    const int report_step_idx = this->ebos_simulator_.episodeIndex();
+    auto ecl_well = this->std_well_.wellEcl();
+    this->well_name_ = ecl_well.name();
+    const GasLiftOpt& glo = schedule.glo(report_step_idx);
+    // NOTE: According to LIFTOPT, item 1:
+    //   "Increment size for lift gas injection rate. Lift gas is
+    //   allocated to individual wells in whole numbers of the increment
+    //   size. If gas lift optimization is no longer required, it can be
+    //   turned off by entering a zero or negative number."
+    // NOTE: This condition was checked in doGasLiftOptimize() in StandardWell
+    //   so it can be assumed that increment_ > 0
+    this->increment_ = glo.gaslift_increment();
+    assert( this->increment_ > 0);
+    // NOTE: The manual (see LIFTOPT, item 2) does not mention
+    //  any default value or restrictions on the economic gradient.
+    // TODO: The value of the gradient would most likely be a positive
+    //  number. Should we warn or fail on a negative value?
+    //  A negative value for the economic gradient would mean that
+    //  the oil production is decreasing with increased liftgas
+    //  injection (which seems strange)
+    this->eco_grad_ = glo.min_eco_gradient();
+    auto& gl_well = glo.well(this->well_name_);
+
+    if(useFixedAlq_(gl_well)) {
+        updateWellStateAlqFixedValue_(gl_well);
+        this->optimize_ = false; // lift gas supply is fixed
+    }
+    else {
+        setAlqMaxRate_(gl_well);
+        this->optimize_ = true;
+    }
+    computeInitialWellRates_();
+
+    if(this->optimize_) {
+        setAlqMinRate_(gl_well);
+        // NOTE: According to item 4 in WLIFTOPT, this value does not
+        //    have to be positive.
+        // TODO: Does it make sense to have a negative value?
+        this->alpha_w_ = gl_well.weight_factor();
+        if (this->alpha_w_ <= 0 ) {
+            displayWarning_("Nonpositive value for alpha_w ignored");
+            this->alpha_w_ = 1.0;
+        }
+
+        // NOTE: According to item 6 in WLIFTOPT:
+        //   "If this value is greater than zero, the incremental gas rate will influence
+        //    the calculation of the incremental gradient and may be used
+        //    to discourage the allocation of lift gas to wells which produce more gas."
+        // TODO: Does this mean that we should ignore this value if it
+        //   is negative?
+        this->alpha_g_ = gl_well.inc_weight_factor();
+
+        const auto& pu = std_well_.phaseUsage();
+        this->oil_pos_ = pu.phase_pos[Oil];
+        this->gas_pos_ = pu.phase_pos[Gas];
+        this->water_pos_ = pu.phase_pos[Water];
+        this->new_alq_ = this->orig_alq_;
+        // TODO: adhoc value.. Should we keep max_iterations_ as a safety measure
+        //   or does it not make sense to have it?
+        this->max_iterations_ = 1000;
+    }
+}
+
+/****************************************
+ * Methods in alphabetical order
+ ****************************************/
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+computeInitialWellRates_()
+{
+    // get the alq value used for this well for the previous time step, or
+    //   if gas lift optimization has not been applied to this well yet, the
+    //   default value is used.
+    this->orig_alq_ = this->well_state_.getALQ(this->well_name_);
+    // NOTE: compute initial rates with current ALQ
+    this->std_well_.computeWellRatesWithThpAlqProd(
+        this->ebos_simulator_, this->summary_state_, this->deferred_logger_,
+        this->potentials_, this->orig_alq_);
+
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+computeWellRates_(double bhp, std::vector<double> &potentials)
+{
+    this->std_well_.computeWellRatesWithBhp(
+        this->ebos_simulator_, bhp, potentials, this->deferred_logger_);
+
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+debugShowIterationInfo_(OptimizeState &state, double alq)
+{
+    const std::string msg = fmt::format("iteration {}, ALQ = {}", state.it, alq);
+    this->displayDebugMessage_(msg);
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+debugShowStartIteration_(double alq, bool increase)
+{
+    const std::string msg =
+        fmt::format("starting {} iteration, ALQ = {}, oilrate = {}",
+            (increase ? "increase" : "decrease"),
+            alq,
+            -this->potentials_[this->oil_pos_]);
+    this->displayDebugMessage_(msg);
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+displayDebugMessage_(const std::string &msg)
+{
+
+    const std::string message = fmt::format(
+        "  GLIFT (DEBUG) : Well {} : {}", this->well_name_, msg);
+    this->deferred_logger_.info(message);
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+displayWarning_(std::string msg)
+{
+    const std::string message = fmt::format(
+        "GAS LIFT OPTIMIZATION, WELL {} : {}", this->well_name_, msg);
+    this->deferred_logger_.warning("WARNING", message);
+}
+
+// TODO: what if the gas_rate_target_ has been defaulted
+//   (i.e. value == 0, meaning: "No limit") but the
+//   oil_rate_target_ has not been defaulted ?
+//   If the new_oil_rate exceeds the oil_rate_target_ it is cut back,
+//   but the same cut-back will not happen for the new_gas_rate
+//   Seems like an inconsistency, since alq should in this
+//   case also be adjusted (to the smaller value that would
+//   give oil target rate) but then the gas rate would also be smaller?
+//   The effect of not reducing the gas rate (if it should be
+//   reduced?) is that a too large value is used in the
+//   computation of the economic gradient making the gradient
+//   smaller than it should be since the term appears in the denominator.
+template<typename TypeTag>
+double
+Opm::GasLiftRuntime<TypeTag>::
+getGasRateWithLimit_(std::vector<double> &potentials)
+{
+    auto new_rate = -potentials[this->gas_pos_];
+    if (this->controls_.hasControl(Well::ProducerCMode::GRAT)) {
+        auto target = this->controls_.gas_rate;
+        if (new_rate > target)
+            new_rate = target;
+    }
+    return new_rate;
+}
+
+
+// NOTE: If the computed oil rate is larger than the target
+//   rate of the well, we reduce it to the target rate. This
+//   will make the economic gradient smaller than it would be
+//   if we did not reduce the rate, and it is less
+//   likely that the current gas lift increment will be
+//   accepted.
+// TODO: If it still is accepted, we should ideally reduce the alq
+//  also since we also reduced the rate. This might involve
+//   some sort of iteration though..
+template<typename TypeTag>
+double
+Opm::GasLiftRuntime<TypeTag>::
+getOilRateWithLimit_(std::vector<double> &potentials)
+{
+    auto new_rate = -potentials[this->oil_pos_];
+    if (this->controls_.hasControl(Well::ProducerCMode::ORAT)) {
+        auto target = this->controls_.oil_rate;
+        if (new_rate > target)
+            new_rate = target;
+    }
+    return new_rate;
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+logSuccess_()
+{
+
+    const std::string message = fmt::format(
+         "GLIFT, WELL {} {} ALQ from {} to {}",
+         this->well_name_,
+         ((this->new_alq_ > this->orig_alq_) ? "increased" : "decreased"),
+         this->orig_alq_,
+         this->new_alq_);
+    this->deferred_logger_.info(message);
+}
+
+/* - At this point we know that this is a production well, and that its current
+ * control mode is THP.
+ *
+ * - We would like to check if it is possible to
+ *   1) increase the oil production by adding lift gas injection to the
+ *   well, or if that is not possible, if we 2) should reduce the amount
+ *   of lift gas injected due to a too small gain in oil production
+ *   (with the current lift gas injection rate)
+ * - For 1) above, we should not add lift gas if it would cause an oil
+ *   rate target to be exceeded, and for 2) we should not reduce the
+ *   amount of liftgas injected below the minimum lift gas injection
+ *   rate.
+ *
+ *   NOTE: If reducing or adding lift-gas further would cause
+ *     one of the well targets like ORAT, WRAT, GRAT, LRAT, CRAT, RESV, BHP,
+ *     to become violated we should stop the lift gas optimization
+ *     loop.. and then updateWellControls() will later (hopefully) switch the well's
+ *     control mode from THP to the mode of the violated target.
+ *
+ * - Lift gas is added if it is economical viable depending on
+ * the ratio of oil gained compared to the amount of liftgas added.
+ *
+ * - Lift gas supply may be limited.
+ *
+ * - The current value of liftgas for the well is stored in the WellState object.
+ *
+ * - It is assumed that the oil production rate is concave function F
+ *   of the amount of lift gas, such that it increases initially due to the
+ *   reduced density of the mixture in the tubing. However, as the
+ *   lift gas supply is increased further, friction pressure losses in the
+ *   tubing become more important, and the production rate peaks and
+ *   then starts to decrease.
+ *   Since lift gas injection has a price, e.g. compression costs can
+ *   be expressed as a cost per unit rate of lift gas injection,
+ *   it must be balanced against the value of the extra amount of
+ *   oil produced. Thus there is a  "minimum economic gradient" of oil
+ *   production rate versus lift gas injection rate, at which the
+ *   value of the extra amount of oil produced by a small increase in
+ *   the lift gas injection rate is equal to the cost of supplying the
+ *   extra amount of lift gas. The optimum lift gas injection rate is then somewhat
+ *   lower than the peak value.
+ *
+ *   Based on this assumption, we know that if the gradient (derivative) of F is
+ *   positive, but greater than the economic gradient (assuming the
+ *   economic gradient is positive), we should add
+ *   lift gas. On the other hand, if the gradient of F is negative or
+ *   if it is positive but smaller than the economic gradient, the amount
+ *   of lift gas injected should be decreased.
+ *
+ */
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+runOptimize()
+{
+    if (this->optimize_) {
+        if (!tryIncreaseLiftGas_()) {
+            if (!tryDecreaseLiftGas_()) {
+                return;
+            }
+        }
+        logSuccess_();
+        this->well_state_.setALQ(this->well_name_, this->new_alq_);
+    }
+    // NOTE: In addition to the new ALQ value, we also implicitly
+    //    return this->potentials_
+}
+
+// INPUT:
+// - increase (boolean) :
+//   - true  : try increase the lift gas supply,
+//   - false : try decrease lift gas supply.
+//
+// OUTPUT:
+// - return value: success (true/false)
+// - potentials_ : updated well potentials if success
+// - new_alq_    : updated alq if success
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::
+runOptimizeLoop_(bool increase)
+{
+    auto cur_potentials = this->potentials_;  // make copy, since we may fail..
+    auto oil_rate = -cur_potentials[this->oil_pos_];
+    auto gas_rate = -cur_potentials[this->gas_pos_];
+    bool success = false;  // did we succeed to increase alq?
+    auto cur_alq = this->orig_alq_;
+    auto alq = cur_alq;
+    OptimizeState state {*this, increase};
+    if (this->debug) debugShowStartIteration_(alq, increase);
+    while (!state.stop_iteration && (++state.it <= this->max_iterations_)) {
+        if (state.checkWellRatesViolated(cur_potentials)) break;
+        if (state.checkAlqOutsideLimits(alq, oil_rate)) break;
+        alq = state.addOrSubtractAlqIncrement(alq);
+        if(this->debug) debugShowIterationInfo_(state, alq);
+        if (!state.computeBhpAtThpLimit(alq)) break;
+        // NOTE: if BHP is below limit, we set state.stop_iteration = true
+        auto bhp = state.getBhpWithLimit();
+        computeWellRates_(bhp, cur_potentials);
+        auto new_oil_rate = getOilRateWithLimit_(cur_potentials);
+        auto new_gas_rate = getGasRateWithLimit_(cur_potentials);
+        auto gradient = state.calcGradient(
+            oil_rate, new_oil_rate, gas_rate, new_gas_rate);
+        if (state.checkEcoGradient(gradient)) {
+            if (state.it == 1) {
+                break;
+            }
+            else {
+                state.stop_iteration = true;
+            }
+        }
+        cur_alq = alq;
+        success = true;
+        oil_rate = new_oil_rate;
+        gas_rate = new_gas_rate;
+    }
+    if (state.it > this->max_iterations_) {
+        warnMaxIterationsExceeded_();
+    }
+    if (success) {
+        this->potentials_ = cur_potentials;
+        this->new_alq_ = cur_alq;
+    }
+    return success;
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+setAlqMaxRate_(const GasLiftOpt::Well &well)
+{
+    auto& max_alq_optional = well.max_rate();
+    if (max_alq_optional) {
+        // NOTE: To prevent extrapolation of the VFP tables, any value
+        // entered here must not exceed the largest ALQ value in the well's VFP table.
+        this->max_alq_ = *max_alq_optional;
+    }
+    else { // i.e. WLIFTOPT, item 3 has been defaulted
+        // According to the Eclipse manual for WLIFTOPT, item 3:
+        //   The default value should be set to the largest ALQ
+        //   value in the well's VFP table
+        const auto& table = *(std_well_.vfp_properties_->getProd()->getTable(
+                this->controls_.vfp_table_number));
+        const auto& alq_values = table.getALQAxis();
+        // Assume the alq_values are sorted in ascending order, so
+        // the last item should be the largest value:
+        this->max_alq_ = alq_values.back();
+    }
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+setAlqMinRate_(const GasLiftOpt::Well &well)
+{
+    // NOTE:  According to WLIFTOPT item 5 :
+    //   if min_rate() is negative, it means: allocate at least enough lift gas
+    //   to enable the well to flow
+    // NOTE: "to enable the well to flow" : How to interpret this?
+    //   We choose to interpret it to mean a positive oil rate as returned from
+    //
+    //    computeWellRates_(bhp, cur_potentials);
+    //
+    //   So even if the well is producing gas, if the oil rate is zero
+    //   we say that the "well is not flowing".
+    //
+    //   Note that if WECON item 2 is set, the well can be shut off
+    //   before the flow rate reaches zero. Also,
+    //   if bhp drops below the bhp lower limit, the well might switch to bhp
+    //   control before the oil rate becomes zero.
+
+    this->min_alq_ = well.min_rate();
+    if (this->min_alq_ > 0) {
+        if (this->min_alq_ >= this->max_alq_) {
+            // NOTE: We reset the value to a negative value.
+            //   negative value means: Allocate at least enough lift gas
+            //   to allow the well to flow.
+            // TODO: Consider other options for resetting the value..
+            this->min_alq_ = -1;
+            displayWarning_("Minimum ALQ value is larger than maximum ALQ value!"
+                " Resetting value.");
+        }
+    }
+}
+
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::
+tryDecreaseLiftGas_()
+{
+    return runOptimizeLoop_(/*increase=*/ false);
+}
+
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::
+tryIncreaseLiftGas_()
+{
+    return runOptimizeLoop_(/*increase=*/ true);
+}
+
+
+// Called when we should use a fixed ALQ value
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+updateWellStateAlqFixedValue_(const GasLiftOpt::Well &well)
+{
+    auto& max_alq_optional = well.max_rate();
+    if (max_alq_optional) {
+        // According to WLIFTOPT, item 3:
+        // If item 2 is NO, then item 3 is regarded as the fixed
+        // lift gas injection rate for the well.
+        auto new_alq = *max_alq_optional;
+        this->well_state_.setALQ(this->well_name_, new_alq);
+    }
+    // else {
+    //    // If item 3 is defaulted, the lift gas rate remains
+    //    // unchanged at its current value.
+    //}
+
+}
+
+// Determine if we should use a fixed ALQ value.
+//
+// From the manual for WLIFTOPT, item 2:
+//   Is the well's lift gas injection rate to be calculated by the
+//   optimization facility?
+// - YES : The well's lift gas injection rate is calculated by the
+//   optimization facility.
+// - NO  : The well's lift gas injection rate remains fixed at a
+//   value that can be set either in Item 3 of this keyword, or in
+//   Item 12 of keyword WCONPROD, or with keyword WELTARG.
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::
+useFixedAlq_(const GasLiftOpt::Well &well)
+{
+    auto wliftopt_item2 = well.use_glo();
+    if (wliftopt_item2) {
+        return false;
+    }
+    else {
+        //  auto& max_alq_optional = well.max_rate();
+        //  if (max_alq_optional) {
+               // According to WLIFTOPT, item 3:
+               // If item 2 is NO, then item 3 is regarded as the fixed
+               // lift gas injection rate for the well.
+        //  }
+        //  else {
+              // If item 3 is defaulted, the lift gas rate remains
+              // unchanged at its current value.
+        //  }
+        return true;
+    }
+}
+
+template<typename TypeTag>
+void
+Opm::GasLiftRuntime<TypeTag>::
+warnMaxIterationsExceeded_()
+{
+    std::ostringstream ss;
+    ss << "Max iterations (" << this->max_iterations_ << ") exceeded in "
+       << "gas lift optimization for well " << this->well_name_;
+    deferred_logger_.warning("MAX_ITERATIONS_EXCEEDED", ss.str());
+}
+
+/****************************************
+ * Methods declared in OptimizeState
+ ****************************************/
+
+template<typename TypeTag>
+double
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+addOrSubtractAlqIncrement(double alq)
+{
+    if (this->increase) {
+        alq += this->parent.increment_;
+        // NOTE: if max_alq_ was defaulted in WCONPROD, item 3, it has
+        //   already been set to the largest value in the VFP table in
+        //   the contructor of GasLiftRuntime
+        if (alq > this->parent.max_alq_) alq = this->parent.max_alq_;
+    }
+    else {
+        alq -= this->parent.increment_;
+        if (this->parent.min_alq_ > 0) {
+            if (alq < this->parent.min_alq_) alq = this->parent.min_alq_;
+        }
+    }
+    return alq;
+}
+
+template<typename TypeTag>
+double
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+calcGradient(double oil_rate, double new_oil_rate, double gas_rate, double new_gas_rate)
+{
+    auto dqo = new_oil_rate - oil_rate;
+    auto dqg = new_gas_rate - gas_rate;
+    // TODO: Should we do any error checks on the calculation of the
+    //   gradient?
+    // NOTE: The eclipse techincal description (chapter 25) says:
+    //   "The gas rate term in the denominator is subject to the
+    //   constraint alpha_g_ * dqg >= 0.0"
+    auto gradient = (this->parent.alpha_w_ * dqo) /
+        (this->parent.increment_ + this->parent.alpha_g_*dqg);
+    return gradient;
+}
+
+// NOTE:  According to WLIFTOPT item 5 :
+//   if min_rate() is negative, it means: allocate at least enough lift gas
+//   to enable the well to flow
+//  We will interpret this as (see discussion above GasLiftRuntime()
+//   in this file): Allocate at least the amount of lift gas needed to
+//   get a positive oil production rate.
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+checkAlqOutsideLimits(double alq, double oil_rate)
+{
+    std::ostringstream ss;
+    bool result = false;
+
+    if (this->increase) {
+        if (alq >= this->parent.max_alq_) {
+            ss << "ALQ >= " << this->parent.max_alq_ << " (max limit), "
+               << "stopping iteration";
+            result = true;
+        }
+        else {
+            // NOTE: A negative min_alq_ means: allocate at least enough lift gas
+            //  to enable the well to flow, see WCONPROD item 5.
+            if (this->parent.min_alq_ < 0) {
+                result = false;
+            }
+            else {
+                // NOTE: checking for a lower limit should not be necessary
+                // when increasing alq.. so this is just to catch an
+                // illegal state at an early point.
+                if (alq < this->parent.min_alq_ ) {
+                    warn_("unexpected: alq below lower limit when trying to "
+                        "increase lift gas. aborting iteration.");
+                    result = true;
+                }
+                else {
+                    result = false;
+                }
+            }
+        }
+    }
+    else { // we are decreasing lift gas
+        // NOTE: A negative min_alq_ means: allocate at least enough lift gas
+        //  to enable the well to flow, see WCONPROD item 5.
+        if (this->parent.min_alq_ < 0) {
+            // If the oil rate is already zero or negative (non-flowing well)
+            // we assume we will not be able to increase it by decreasing the lift gas
+            if ( oil_rate <= 0 ) {
+                ss << "Oil rate ( " << oil_rate << " ) <= 0 when decreasing lift gas. "
+                   << "We will not be able to make this well flowing by decreasing "
+                   << "lift gas, stopping iteration.";
+                result = true;
+            }
+            else {
+                result = false;
+            }
+        }
+        else {
+            if (alq <= this->parent.min_alq_ ) {
+                ss << "ALQ <= " << this->parent.min_alq_ << " (min limit), "
+                   << "stopping iteration";
+                result = true;
+            }
+            else {
+                // 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_) {
+                    warn_( "unexpected: alq above upper limit when trying to "
+                        "decrease lift gas. aborting iteration.");
+                    result = true;
+                }
+                else {
+                    result = false;
+                }
+            }
+        }
+    }
+    if (this->parent.debug) this->parent.displayDebugMessage_(ss.str());
+    return result;
+}
+
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+checkEcoGradient(double gradient)
+{
+    std::ostringstream ss;
+    bool result = false;
+
+    if (this->parent.debug) {
+        ss << "checking gradient: " << gradient;
+    }
+    if (this->increase) {
+        if (this->parent.debug) ss << " <= " << this->parent.eco_grad_ << " --> ";
+        if (gradient <= this->parent.eco_grad_) {
+            if (this->parent.debug) ss << "true";
+            result = true;
+        }
+        else {
+            if (this->parent.debug) ss << "false";
+        }
+    }
+    else {  // decreasing lift gas
+        if (this->parent.debug) ss << " >= " << this->parent.eco_grad_ << " --> ";
+        if (gradient >= this->parent.eco_grad_) {
+            if (this->parent.debug) ss << "true";
+            result = true;
+        }
+        else {
+            if (this->parent.debug) ss << "false";
+        }
+    }
+    if (this->parent.debug) this->parent.displayDebugMessage_(ss.str());
+    return result;
+}
+
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+checkRate(double rate, double limit, const std::string rate_str)
+{
+
+    if (limit < rate ) {
+        if (this->parent.debug) {
+            const std::string msg = fmt::format(
+                "iteration {} : rate {} exceeds target rate {}. Stopping iteration",
+                this->it, rate_str, rate, limit);
+            this->parent.displayDebugMessage_(msg);
+        }
+        return true;
+    }
+    return false;
+}
+
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+checkWellRatesViolated(std::vector<double> &potentials)
+{
+    if (this->parent.controls_.hasControl(Well::ProducerCMode::ORAT)) {
+        auto oil_rate = -potentials[this->parent.oil_pos_];
+        if (this->checkRate(oil_rate, this->parent.controls_.oil_rate, "oil"))
+            return true;
+    }
+    if (this->parent.controls_.hasControl(Well::ProducerCMode::WRAT)) {
+        auto water_rate = -potentials[this->parent.water_pos_];
+        if (this->checkRate(water_rate, this->parent.controls_.water_rate, "water"))
+            return true;
+    }
+    if (this->parent.controls_.hasControl(Well::ProducerCMode::GRAT)) {
+        auto gas_rate = -potentials[this->parent.gas_pos_];
+        if (this->checkRate(gas_rate, this->parent.controls_.gas_rate, "gas"))
+            return true;
+    }
+    if (this->parent.controls_.hasControl(Well::ProducerCMode::LRAT)) {
+        auto oil_rate = -potentials[this->parent.oil_pos_];
+        auto water_rate = -potentials[this->parent.water_pos_];
+        auto liq_rate = oil_rate + water_rate;
+        if (this->checkRate(liq_rate, this->parent.controls_.liquid_rate, "liquid"))
+            return true;
+    }
+    // TODO: Also check RESV, see checkIndividualContraints() in
+    //   WellInterface_impl.hpp
+    // TODO: Check group contraints?
+
+    return false;
+}
+
+
+
+template<typename TypeTag>
+bool
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+computeBhpAtThpLimit(double alq)
+{
+    auto bhp_at_thp_limit = this->parent.std_well_.computeBhpAtThpLimitProdWithAlq(
+        this->parent.ebos_simulator_,
+        this->parent.summary_state_,
+        this->parent.deferred_logger_,
+        alq);
+    if (!bhp_at_thp_limit) {
+        const std::string msg = fmt::format(
+          "Failed in getting converged bhp potential for well {}",
+          this->parent.well_name_);
+        this->parent.deferred_logger_.warning(
+            "FAILURE_GETTING_CONVERGED_POTENTIAL", msg);
+        return false;
+    }
+    this->bhp = *bhp_at_thp_limit;
+    return true;
+}
+
+//  NOTE: When calculating the gradient, determine what the well would produce if
+//  the lift gas injection rate were increased by one increment. The
+//  production rates are adjusted if necessary to obey
+//  any rate or BHP limits that the well may be subject to. From this
+//  information, calculate the well's "weighted incremental
+//  gradient"
+//
+// TODO: What does it mean to "adjust the production rates" given a
+//   BHP limit?
+//
+template<typename TypeTag>
+double
+Opm::GasLiftRuntime<TypeTag>::OptimizeState::
+getBhpWithLimit()
+{
+    auto bhp_update = this->bhp;
+    if (this->parent.controls_.hasControl(Well::ProducerCMode::BHP)) {
+        auto limit = this->parent.controls_.bhp_limit;
+        // TODO: is it possible that bhp falls below the limit when
+        // adding lift gas? I.e. if this->increase == true..
+        if (this->bhp < limit) {
+            // TODO: we keep the current alq, but it should probably
+            // be adjusted since we changed computed bhp. But how?
+            bhp_update = limit;
+            // Stop iteration, but first check the economic gradient
+            //   with the bhp_update. If the gradient looks OK (see the
+            //   main optimize loop) we keep the current ALQ value.
+            this->stop_iteration = true;
+        }
+    }
+    return bhp_update;
+}

opm/simulators/wells/MultisegmentWell.hpp

@@ -115,6 +115,14 @@ namespace Opm
 
         virtual void initPrimaryVariablesEvaluation() const override;
 
+        virtual void maybeDoGasLiftOptimization (
+            const WellState&,
+            const Simulator&,
+            DeferredLogger&
+        ) const {
+            // Not implemented yet
+        }
+
         virtual void assembleWellEq(const Simulator& ebosSimulator,
                                     const std::vector<Scalar>& B_avg,
                                     const double dt,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -2169,8 +2169,8 @@ namespace Opm
             const EvalWell inlet_mass_rate = segment_mass_rates_[inlet];
             accelerationPressureLoss -= mswellhelpers::velocityHead(std::max(inlet_area, area), inlet_mass_rate, inlet_density);
         }
-        
-        // We change the sign of the accelerationPressureLoss for injectors. 
+
+        // We change the sign of the accelerationPressureLoss for injectors.
         // Is this correct? Testing indicates that this is what the reference simulator does
         const double sign = mass_rate < 0. ? 1.0 : - 1.0;
         accelerationPressureLoss *= sign;
@@ -2539,7 +2539,7 @@ namespace Opm
                 const EvalWell segment_surface_volume = getSegmentSurfaceVolume(ebosSimulator, seg);
 
                 // Add a regularization_factor to increase the accumulation term
-                // This will make the system less stiff and help convergence for 
+                // This will make the system less stiff and help convergence for
                 // difficult cases
                 const Scalar regularization_factor =  param_.regularization_factor_ms_wells_;
                 // for each component

opm/simulators/wells/StandardWell.hpp

@@ -29,6 +29,7 @@
 
 #include <opm/simulators/wells/RateConverter.hpp>
 #include <opm/simulators/wells/WellInterface.hpp>
+#include <opm/simulators/wells/GasLiftRuntime.hpp>
 
 #include <opm/models/blackoil/blackoilpolymermodules.hh>
 #include <opm/models/blackoil/blackoilsolventmodules.hh>
@@ -42,6 +43,7 @@
 #include <dune/common/dynmatrix.hh>
 
 #include <optional>
+#include <fmt/format.h>
 
 namespace Opm
 {
@@ -66,6 +68,7 @@ namespace Opm
         using typename Base::RateConverterType;
         using typename Base::SparseMatrixAdapter;
         using typename Base::FluidState;
+        using GasLiftHandler = Opm::GasLiftRuntime<TypeTag>;
 
         using Base::numEq;
 
@@ -233,11 +236,70 @@ namespace Opm
             return param_.matrix_add_well_contributions_;
         }
 
+        bool doGasLiftOptimize(
+            const WellState& well_state,
+            const Simulator& ebosSimulator,
+            DeferredLogger& deferred_logger
+        ) const;
+
+        virtual void maybeDoGasLiftOptimization (
+            const WellState& well_state,
+            const Simulator& ebosSimulator,
+            DeferredLogger& deferred_logger
+        ) const;
+
+        bool checkGliftNewtonIterationIdxOk(
+            const Simulator& ebosSimulator,
+            DeferredLogger& deferred_logger
+        ) const;
+
+        void gliftDebug(
+            const std::string &msg,
+            Opm::DeferredLogger& deferred_logger) const;
+
+        void gasLiftOptimizeProduction(
+            const Simulator& ebosSimulator,
+            const SummaryState& summaryState,
+            DeferredLogger& deferredLogger,
+            std::vector<double>& potentials,
+            const WellState& well_state);
+
+        /* returns BHP */
+        double computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
+                               const SummaryState &summary_state,
+                               DeferredLogger &deferred_logger,
+                               std::vector<double> &potentials,
+                               double alq) const;
+
+        void computeWellRatesWithThpAlqProd(
+            const Simulator &ebos_simulator,
+            const SummaryState &summary_state,
+            DeferredLogger &deferred_logger,
+            std::vector<double> &potentials,
+            double alq) const;
+
+        // NOTE: Cannot be protected since it is used by GasLiftRuntime
+        std::optional<double> computeBhpAtThpLimitProdWithAlq(
+            const Simulator& ebos_simulator,
+            const SummaryState& summary_state,
+            DeferredLogger& deferred_logger,
+            double alq_value) const;
+
+        // NOTE: Cannot be protected since it is used by GasLiftRuntime
+        void computeWellRatesWithBhp(
+            const Simulator& ebosSimulator,
+            const double& bhp,
+            std::vector<double>& well_flux,
+            Opm::DeferredLogger& deferred_logger) const;
+
+        // NOTE: These cannot be protected since they are used by GasLiftRuntime
+        using Base::phaseUsage;
+        using Base::vfp_properties_;
+
     protected:
 
         // protected functions from the Base class
         using Base::getAllowCrossFlow;
-        using Base::phaseUsage;
         using Base::flowPhaseToEbosCompIdx;
         using Base::ebosCompIdxToFlowCompIdx;
         using Base::wsalt;
@@ -251,7 +313,6 @@ namespace Opm
         // protected member variables from the Base class
         using Base::current_step_;
         using Base::well_ecl_;
-        using Base::vfp_properties_;
         using Base::gravity_;
         using Base::param_;
         using Base::well_efficiency_factor_;
@@ -313,6 +374,8 @@ namespace Opm
         mutable std::vector<double> ipr_b_;
 
         bool changed_to_stopped_this_step_ = false;
+        // Enable GLIFT debug mode. This will enable output of logging messages.
+        bool glift_debug = false;
 
         const EvalWell& getBhp() const;
 
@@ -380,25 +443,22 @@ namespace Opm
                              double& perf_vap_oil_rate,
                              Opm::DeferredLogger& deferred_logger) const;
 
-        void computeWellRatesWithBhp(const Simulator& ebosSimulator,
-                                             const double& bhp,
-                                             std::vector<double>& well_flux,
-                                             Opm::DeferredLogger& deferred_logger) const;
-
         void computeWellRatesWithBhpPotential(const Simulator& ebosSimulator,
                                               const std::vector<Scalar>& B_avg,
                                               const double& bhp,
                                               std::vector<double>& well_flux,
                                               Opm::DeferredLogger& deferred_logger);
 
-        std::vector<double> computeWellPotentialWithTHP(const Simulator& ebosSimulator,
-                                                        Opm::DeferredLogger& deferred_logger) const;
+        std::vector<double> computeWellPotentialWithTHP(
+            const Simulator& ebosSimulator,
+            Opm::DeferredLogger& deferred_logger,
+            const WellState &well_state) const;
 
         template <class ValueType>
-        ValueType calculateBhpFromThp(const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger) const;
+        ValueType calculateBhpFromThp(const WellState& well_state, const std::vector<ValueType>& rates, const Well& well, const SummaryState& summaryState, Opm::DeferredLogger& deferred_logger) const;
 
 
-        double calculateThpFromBhp(const std::vector<double>& rates, const double bhp, Opm::DeferredLogger& deferred_logger) const;
+        double calculateThpFromBhp(const WellState &well_state, const std::vector<double>& rates, const double bhp, Opm::DeferredLogger& deferred_logger) const;
 
         // get the mobility for specific perforation
         void getMobility(const Simulator& ebosSimulator,
@@ -422,6 +482,8 @@ namespace Opm
 
         void updateThp(WellState& well_state, Opm::DeferredLogger& deferred_logger) const;
 
+        double getALQ(const WellState& well_state) const;
+
         void assembleControlEq(const WellState& well_state,
                                const Opm::Schedule& schedule,
                                const SummaryState& summaryState,
@@ -453,10 +515,10 @@ namespace Opm
                                    Opm::DeferredLogger& deferred_logger);
 
         // check whether the well is operable under BHP limit with current reservoir condition
-        void checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
+        void checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
 
         // check whether the well is operable under THP limit with current reservoir condition
-        void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger);
+        void checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger);
 
         // for a well, when all drawdown are in the wrong direction, then this well will not
         // be able to produce/inject .
@@ -537,7 +599,8 @@ namespace Opm
                                         DeferredLogger& deferred_logger);
 
 
-        std::optional<double> computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
+        std::optional<double> computeBhpAtThpLimitProd(const WellState& well_state,
+                                                       const Simulator& ebos_simulator,
                                                        const SummaryState& summary_state,
                                                        DeferredLogger& deferred_logger) const;
 

opm/simulators/wells/StandardWell_impl.hpp

@@ -827,7 +827,7 @@ namespace Opm
             // Setup function for evaluation of BHP from THP (used only if needed).
             auto bhp_from_thp = [&]() {
                 const auto rates = getRates();
-                return calculateBhpFromThp(rates, well, summaryState, deferred_logger);
+                return calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
             };
             // Call generic implementation.
             const auto& inj_controls = well.injectionControls(summaryState);
@@ -837,7 +837,7 @@ namespace Opm
             const auto rates = getRates();
             // Setup function for evaluation of BHP from THP (used only if needed).
             auto bhp_from_thp = [&]() {
-                return calculateBhpFromThp(rates, well, summaryState, deferred_logger);
+                 return calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
             };
             // Call generic implementation.
             const auto& prod_controls = well.productionControls(summaryState);
@@ -1235,7 +1235,7 @@ namespace Opm
 
         const double bhp = well_state.bhp()[index_of_well_];
 
-        well_state.thp()[index_of_well_] = calculateThpFromBhp(rates, bhp, deferred_logger);
+        well_state.thp()[index_of_well_] = calculateThpFromBhp(well_state, rates, bhp, deferred_logger);
 
     }
 
@@ -1316,7 +1316,7 @@ namespace Opm
                 for (int p = 0; p<np; ++p) {
                     rates[p] = well_state.wellRates()[well_index*np + p];
                 }
-                double bhp = calculateBhpFromThp(rates, well, summaryState, deferred_logger);
+                double bhp = calculateBhpFromThp(well_state, rates, well, summaryState, deferred_logger);
                 well_state.bhp()[well_index] = bhp;
                 break;
             }
@@ -1442,7 +1442,10 @@ namespace Opm
             case Well::ProducerCMode::THP:
             {
                 well_state.thp()[well_index] = controls.thp_limit;
-                auto bhp = computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger);
+                gliftDebug(
+                    "computing BHP from THP to update well state",
+                    deferred_logger);
+                auto bhp = computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger);
                 if (bhp) {
                     well_state.bhp()[well_index] = *bhp;
                 } else {
@@ -1628,7 +1631,7 @@ namespace Opm
     void
     StandardWell<TypeTag>::
     updateWellOperability(const Simulator& ebos_simulator,
-                          const WellState& /* well_state */,
+                          const WellState& well_state,
                           Opm::DeferredLogger& deferred_logger)
     {
         this->operability_status_.reset();
@@ -1636,13 +1639,13 @@ namespace Opm
         updateIPR(ebos_simulator, deferred_logger);
 
         // checking the BHP limit related
-        checkOperabilityUnderBHPLimitProducer(ebos_simulator, deferred_logger);
+        checkOperabilityUnderBHPLimitProducer(well_state, ebos_simulator, deferred_logger);
 
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
 
         // checking whether the well can operate under the THP constraints.
         if (this->wellHasTHPConstraints(summaryState)) {
-            checkOperabilityUnderTHPLimitProducer(ebos_simulator, deferred_logger);
+            checkOperabilityUnderTHPLimitProducer(ebos_simulator, well_state, deferred_logger);
         }
     }
 
@@ -1653,7 +1656,7 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    checkOperabilityUnderBHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger)
+    checkOperabilityUnderBHPLimitProducer(const WellState& well_state, const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger)
     {
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
         const double bhp_limit = mostStrictBhpFromBhpLimits(summaryState);
@@ -1680,7 +1683,7 @@ namespace Opm
                 std::vector<double> well_rates_bhp_limit;
                 computeWellRatesWithBhp(ebos_simulator, bhp_limit, well_rates_bhp_limit, deferred_logger);
 
-                const double thp = calculateThpFromBhp(well_rates_bhp_limit, bhp_limit, deferred_logger);
+                const double thp = calculateThpFromBhp(well_state, well_rates_bhp_limit, bhp_limit, deferred_logger);
                 const double thp_limit = this->getTHPConstraint(summaryState);
 
                 if (thp < thp_limit) {
@@ -1707,10 +1710,10 @@ namespace Opm
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::
-    checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, Opm::DeferredLogger& deferred_logger)
+    checkOperabilityUnderTHPLimitProducer(const Simulator& ebos_simulator, const WellState& well_state, Opm::DeferredLogger& deferred_logger)
     {
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
-        const auto obtain_bhp = computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger);
+        const auto obtain_bhp = computeBhpAtThpLimitProd(well_state, ebos_simulator, summaryState, deferred_logger);
 
         if (obtain_bhp) {
             this->operability_status_.can_obtain_bhp_with_thp_limit = true;
@@ -2507,48 +2510,214 @@ namespace Opm
     std::vector<double>
     StandardWell<TypeTag>::
     computeWellPotentialWithTHP(const Simulator& ebos_simulator,
-                                Opm::DeferredLogger& deferred_logger) const
+                               Opm::DeferredLogger& deferred_logger,
+                               const WellState &well_state) const
     {
         std::vector<double> potentials(number_of_phases_, 0.0);
         const auto& summary_state = ebos_simulator.vanguard().summaryState();
 
         const auto& well = well_ecl_;
         if (well.isInjector()){
+            const auto& controls = well_ecl_.injectionControls(summary_state);
             auto bhp_at_thp_limit = computeBhpAtThpLimitInj(ebos_simulator, summary_state, deferred_logger);
             if (bhp_at_thp_limit) {
-                const auto& controls = well_ecl_.injectionControls(summary_state);
                 const double bhp = std::min(*bhp_at_thp_limit, controls.bhp_limit);
                 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
             } else {
                 deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
                                         "Failed in getting converged thp based potential calculation for well "
                                         + name() + ". Instead the bhp based value is used");
-                const auto& controls = well_ecl_.injectionControls(summary_state);
                 const double bhp = controls.bhp_limit;
                 computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
             }
         } else {
-            auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
-            if (bhp_at_thp_limit) {
-                const auto& controls = well_ecl_.productionControls(summary_state);
-                const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
-                computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
-            } else {
-                deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
-                                        "Failed in getting converged thp based potential calculation for well "
-                                        + name() + ". Instead the bhp based value is used");
-                const auto& controls = well_ecl_.productionControls(summary_state);
-                const double bhp = controls.bhp_limit;
-                computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
-            }
+            computeWellRatesWithThpAlqProd(
+                ebos_simulator, summary_state,
+                deferred_logger, potentials, getALQ(well_state)
+            );
         }
 
         return potentials;
     }
 
+    template<typename TypeTag>
+    double
+    StandardWell<TypeTag>::
+    getALQ(const WellState &well_state) const
+    {
+        const auto& well = well_ecl_;
+        const std::string &well_name = well.name();
+        return well_state.getALQ(well_name);
+    }
+
+    /* At this point we know that the well does not have BHP control mode and
+       that it does have THP constraints, see computeWellPotentials().
+     * TODO: Currently we limit the application of gas lift optimization to wells
+     *    operating under THP control mode, does it make sense to
+     *    extend it to other modes?
+     */
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    doGasLiftOptimize(const WellState &well_state, const Simulator &ebos_simulator,
+                      Opm::DeferredLogger& deferred_logger) const
+    {
+
+        gliftDebug("checking if GLIFT should be done..", deferred_logger);
+        if (!well_state.gliftOptimizationEnabled()) {
+            gliftDebug("Optimization disabled in WellState", deferred_logger);
+            return false;
+        }
+        const int well_index = index_of_well_;
+        const Well::ProducerCMode& control_mode
+            = well_state.currentProductionControls()[well_index];
+        if (control_mode != Well::ProducerCMode::THP ) {
+            gliftDebug("Not THP control", deferred_logger);
+            return false;
+        }
+        if (!checkGliftNewtonIterationIdxOk(ebos_simulator, deferred_logger)) {
+            return false;
+        }
+        const int report_step_idx = ebos_simulator.episodeIndex();
+        const Opm::Schedule& schedule = ebos_simulator.vanguard().schedule();
+        const GasLiftOpt& glo = schedule.glo(report_step_idx);
+        auto increment = glo.gaslift_increment();
+        // NOTE: According to the manual: LIFTOPT, item 1, :
+        //   "Increment size for lift gas injection rate. Lift gas is
+        //   allocated to individual wells in whole numbers of the increment
+        //   size.  If gas lift optimization is no longer required, it can be
+        //   turned off by entering a zero or negative number."
+        if (increment <= 0) {
+            if (this->glift_debug) {
+                const std::string msg = fmt::format(
+                    "Gas Lift switched off in LIFTOPT item 1 due to non-positive "
+                    "value: {}", increment);
+                gliftDebug(msg, deferred_logger);
+            }
+            return false;
+        }
+        else {
+            return true;
+        }
+    }
+
+    template<typename TypeTag>
+    bool
+    StandardWell<TypeTag>::
+    checkGliftNewtonIterationIdxOk(
+            const Simulator& ebos_simulator,
+            DeferredLogger& deferred_logger ) const
+    {
+        const int report_step_idx = ebos_simulator.episodeIndex();
+        const Opm::Schedule& schedule = ebos_simulator.vanguard().schedule();
+        const GasLiftOpt& glo = schedule.glo(report_step_idx);
+        const int iteration_idx = ebos_simulator.model().newtonMethod().numIterations();
+        if (glo.all_newton()) {
+            const int nupcol = schedule.getNupcol(report_step_idx);
+            if (this->glift_debug) {
+                const std::string msg = fmt::format(
+                    "LIFTOPT item4 == YES, it = {}, nupcol = {} -->  GLIFT optimize = {}",
+                    iteration_idx,
+                    nupcol,
+                    ((iteration_idx <= nupcol) ? "TRUE" : "FALSE"));
+                gliftDebug(msg, deferred_logger);
+            }
+            return iteration_idx <= nupcol;
+        }
+        else {
+            if (this->glift_debug) {
+                const std::string msg = fmt::format(
+                    "LIFTOPT item4 == NO, it = {} --> GLIFT optimize = {}",
+                    iteration_idx, ((iteration_idx == 1) ? "TRUE" : "FALSE"));
+                gliftDebug(msg, deferred_logger);
+            }
+            return iteration_idx == 1;
+        }
+    }
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    gliftDebug(
+        const std::string &msg, DeferredLogger& deferred_logger) const
+    {
+        if (this->glift_debug) {
+            const std::string message = fmt::format(
+                "  GLIFT (DEBUG) : Well {} : {}", this->name(), msg);
+            deferred_logger.info(message);
+        }
+    }
+
+    template<typename TypeTag>
+    double
+    StandardWell<TypeTag>::
+    computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
+                               const SummaryState &summary_state,
+                               DeferredLogger &deferred_logger,
+                               std::vector<double> &potentials,
+                               double alq) const
+    {
+        double bhp;
+        auto bhp_at_thp_limit = computeBhpAtThpLimitProdWithAlq(
+                              ebos_simulator, summary_state, deferred_logger, alq);
+        if (bhp_at_thp_limit) {
+            const auto& controls = well_ecl_.productionControls(summary_state);
+            bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
+            computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
+        }
+        else {
+            deferred_logger.warning("FAILURE_GETTING_CONVERGED_POTENTIAL",
+                "Failed in getting converged thp based potential calculation for well "
+                + name() + ". Instead the bhp based value is used");
+            const auto& controls = well_ecl_.productionControls(summary_state);
+            bhp = controls.bhp_limit;
+            computeWellRatesWithBhp(ebos_simulator, bhp, potentials, deferred_logger);
+        }
+        return bhp;
+    }
+
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    computeWellRatesWithThpAlqProd(const Simulator &ebos_simulator,
+                               const SummaryState &summary_state,
+                               DeferredLogger &deferred_logger,
+                               std::vector<double> &potentials,
+                               double alq) const
+    {
+        /*double bhp =*/ computeWellRatesAndBhpWithThpAlqProd(
+                        ebos_simulator, summary_state,
+                        deferred_logger, potentials, alq
+                );
+    }
 
 
-
+    template<typename TypeTag>
+    void
+    StandardWell<TypeTag>::
+    maybeDoGasLiftOptimization(
+                          const WellState& well_state,
+                          const Simulator& ebos_simulator,
+                          Opm::DeferredLogger& deferred_logger) const
+    {
+        const auto& well = well_ecl_;
+        if (well.isProducer()) {
+            const auto& summary_state = ebos_simulator.vanguard().summaryState();
+            const Well::ProducerCMode& current_control
+                = well_state.currentProductionControls()[this->index_of_well_];
+            if ( this->Base::wellHasTHPConstraints(summary_state)
+                && current_control != Well::ProducerCMode::BHP ) {
+                std::vector<double> potentials = well_state.wellPotentials();
+                if (doGasLiftOptimize(well_state, ebos_simulator, deferred_logger)) {
+                    const auto& controls = well.productionControls(summary_state);
+                    GasLiftHandler glift {
+                        *this, ebos_simulator, summary_state,
+                        deferred_logger, potentials, well_state, controls };
+                    glift.runOptimize();
+                }
+            }
+        }
+    }
 
     template<typename TypeTag>
     void
@@ -2581,7 +2750,7 @@ namespace Opm
             well.computeWellRatesWithBhpPotential(ebosSimulator, B_avg, bhp, well_potentials, deferred_logger);
         } else {
             // the well has a THP related constraint
-            well_potentials = well.computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
+            well_potentials = well.computeWellPotentialWithTHP(ebosSimulator, deferred_logger, well_state);
         }
     }
 
@@ -2704,7 +2873,8 @@ namespace Opm
     template<class ValueType>
     ValueType
     StandardWell<TypeTag>::
-    calculateBhpFromThp(const std::vector<ValueType>& rates,
+    calculateBhpFromThp(const WellState &well_state,
+                        const std::vector<ValueType>& rates,
                         const Well& well,
                         const SummaryState& summaryState,
                         Opm::DeferredLogger& deferred_logger) const
@@ -2736,7 +2906,7 @@ namespace Opm
              const auto& controls = well.productionControls(summaryState);
              const double vfp_ref_depth = vfp_properties_->getProd()->getTable(controls.vfp_table_number)->getDatumDepth();
              const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
-             return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit, controls.alq_value) - dp;
+             return vfp_properties_->getProd()->bhp(controls.vfp_table_number, aqua, liquid, vapour, controls.thp_limit, getALQ(well_state)) - dp;
          }
          else {
              OPM_DEFLOG_THROW(std::logic_error, "Expected INJECTOR or PRODUCER well", deferred_logger);
@@ -2751,7 +2921,7 @@ namespace Opm
     template<typename TypeTag>
     double
     StandardWell<TypeTag>::
-    calculateThpFromBhp(const std::vector<double>& rates,
+    calculateThpFromBhp(const WellState &well_state, const std::vector<double>& rates,
                         const double bhp,
                         Opm::DeferredLogger& deferred_logger) const
     {
@@ -2774,7 +2944,7 @@ namespace Opm
         }
         else if (this->isProducer()) {
             const int table_id = well_ecl_.vfp_table_number();
-            const double alq = well_ecl_.alq_value();
+            const double alq = getALQ(well_state);
             const double vfp_ref_depth = vfp_properties_->getProd()->getTable(table_id)->getDatumDepth();
             const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
 
@@ -3402,9 +3572,22 @@ namespace Opm
     template<typename TypeTag>
     std::optional<double>
     StandardWell<TypeTag>::
-    computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
+    computeBhpAtThpLimitProd(const WellState& well_state,
+                             const Simulator& ebos_simulator,
                              const SummaryState& summary_state,
                              DeferredLogger& deferred_logger) const
+    {
+        return computeBhpAtThpLimitProdWithAlq(
+            ebos_simulator, summary_state, deferred_logger, getALQ(well_state));
+    }
+
+    template<typename TypeTag>
+    std::optional<double>
+    StandardWell<TypeTag>::
+    computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
+                             const SummaryState& summary_state,
+                                   DeferredLogger& deferred_logger,
+                                   double alq_value) const
     {
         // Given a VFP function returning bhp as a function of phase
         // rates and thp:
@@ -3451,10 +3634,10 @@ namespace Opm
         const double vfp_ref_depth = table.getDatumDepth();
         const double rho = perf_densities_[0]; // Use the density at the top perforation.
         const double dp = wellhelpers::computeHydrostaticCorrection(ref_depth_, vfp_ref_depth, rho, gravity_);
-        auto fbhp = [this, &controls, dp](const std::vector<double>& rates) {
+        auto fbhp = [this, &controls, dp, alq_value](const std::vector<double>& rates) {
             assert(rates.size() == 3);
             return this->vfp_properties_->getProd()
-            ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit, controls.alq_value) - dp;
+            ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], controls.thp_limit, alq_value) - dp;
         };
 
         // Make the flo() function.

opm/simulators/wells/WellInterface.hpp

@@ -26,6 +26,7 @@
 
 #include <opm/common/OpmLog/OpmLog.hpp>
 #include <opm/common/ErrorMacros.hpp>
+#include <opm/common/Exceptions.hpp>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
 #include <opm/parser/eclipse/EclipseState/Schedule/Well/WellTestState.hpp>
@@ -167,6 +168,12 @@ namespace Opm
                                     Opm::DeferredLogger& deferred_logger
                                     ) = 0;
 
+        virtual void maybeDoGasLiftOptimization (
+            const WellState& well_state,
+            const Simulator& ebosSimulator,
+            DeferredLogger& deferred_logger
+        ) const = 0;
+
         void updateWellTestState(const WellState& well_state,
                                  const double& simulationTime,
                                  const bool& writeMessageToOPMLog,

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -297,6 +297,8 @@ namespace Opm
                 seg_pressdrop_friction_.assign(nw, 0.);
                 seg_pressdrop_acceleration_.assign(nw, 0.);
             }
+            updateWellsDefaultALQ(wells_ecl);
+            do_glift_optimization_ = true;
         }
 
 
@@ -365,7 +367,7 @@ namespace Opm
 
             return it->second;
         }
-        
+
         void setCurrentWellRates(const std::string& wellName, const std::vector<double>& rates ) {
             well_rates[wellName] = rates;
         }
@@ -399,7 +401,7 @@ namespace Opm
         bool hasProductionGroupRates(const std::string& groupName) const {
             return this->production_group_rates.find(groupName) != this->production_group_rates.end();
         }
-        
+
         void setCurrentProductionGroupReductionRates(const std::string& groupName, const std::vector<double>& target ) {
             production_group_reduction_rates[groupName] = target;
         }
@@ -586,6 +588,13 @@ namespace Opm
                     well.rates.set( rt::brine, brineWellRate(w) );
                 }
 
+                if ( well.current_control.isProducer ) {
+                    well.rates.set( rt::alq, getALQ(/*wellName=*/wt.first) );
+                }
+                else {
+                    well.rates.set( rt::alq, 0.0 );
+                }
+
                 well.rates.set( rt::dissolved_gas, this->well_dissolved_gas_rates_[w] );
                 well.rates.set( rt::vaporized_oil, this->well_vaporized_oil_rates_[w] );
 
@@ -1133,6 +1142,41 @@ namespace Opm
             return globalIsProductionGrup_[it->second] != 0;
         }
 
+        void updateALQ( const WellStateFullyImplicitBlackoil &copy ) const
+        {
+            this->current_alq_ = copy.getCurrentALQ();
+        }
+
+        std::map<std::string, double> getCurrentALQ() const
+        {
+            return current_alq_;
+        }
+
+        double getALQ( const std::string& name) const
+        {
+            if (this->current_alq_.count(name) == 0) {
+                this->current_alq_[name] = this->default_alq_[name];
+            }
+            return this->current_alq_[name];
+        }
+
+        void setALQ( const std::string& name, double value) const
+        {
+            this->current_alq_[name] = value;
+        }
+
+        bool gliftOptimizationEnabled() const {
+            return do_glift_optimization_;
+        }
+
+        void disableGliftOptimization() const {
+            do_glift_optimization_ = false;
+        }
+
+        void enableGliftOptimization() const {
+            do_glift_optimization_ = true;
+        }
+
     private:
         std::vector<double> perfphaserates_;
 
@@ -1160,6 +1204,9 @@ namespace Opm
         std::map<std::string, double> injection_group_vrep_rates;
         std::map<std::string, std::vector<double>> injection_group_rein_rates;
         std::map<std::string, double> group_grat_target_from_sales;
+        mutable std::map<std::string, double> current_alq_;
+        mutable std::map<std::string, double> default_alq_;
+        mutable bool do_glift_optimization_;
 
         std::vector<double> perfRateSolvent_;
 
@@ -1289,6 +1336,39 @@ namespace Opm
 
             return this->seg_number_[top_offset + seg_id];
         }
+
+        // If the ALQ has changed since the previous report step,
+        // reset current_alq and update default_alq. ALQ is used for
+        // constant lift gas injection and for gas lift optimization
+        // (THP controlled wells).
+        //
+        // NOTE: If a well is no longer used (e.g. it is shut down)
+        // it is still kept in the maps "default_alq_" and "current_alq_". Since the
+        // number of unused entries should be small (negligible memory
+        // overhead) this is simpler than writing code to delete it.
+        //
+        void updateWellsDefaultALQ( const std::vector<Well>& wells_ecl )
+        {
+            const int nw = wells_ecl.size();
+            for (int i = 0; i<nw; i++) {
+                const Well &well = wells_ecl[i];
+                if (well.isProducer()) {
+                    const std::string &name = well.name();
+                    // NOTE: This is the value set in item 12 of WCONPROD, or with WELTARG
+                    auto alq = well.alq_value();
+                    if (this->default_alq_.count(name) != 0) {
+                        if (this->default_alq_[name] == alq) {
+                            // If the previous value was the same, we leave current_alq_
+                            // as it is.
+                            continue;
+                        }
+                    }
+                    this->default_alq_[name] = alq;
+                    // Reset current ALQ if a new value was given in WCONPROD
+                    this->current_alq_[name] = alq;
+                }
+            }
+        }
     };
 
 } // namespace Opm