Add gaslift optimization support for MSW.

Implements gas lift optimization support for multisegmented wells (MSW).

opm/simulators/wells/BlackoilWellModel.hpp

@@ -105,6 +105,7 @@ namespace Opm {
             using RateVector = GetPropType<TypeTag, Properties::RateVector>;
             using GlobalEqVector = GetPropType<TypeTag, Properties::GlobalEqVector>;
             using SparseMatrixAdapter = GetPropType<TypeTag, Properties::SparseMatrixAdapter>;
+            using GasLiftSingleWell = typename WellInterface<TypeTag>::GasLiftSingleWell;
             using GLiftOptWells = typename BlackoilWellModelGeneric::GLiftOptWells;
             using GLiftProdWells = typename BlackoilWellModelGeneric::GLiftProdWells;
             using GLiftWellStateMap =
@@ -398,6 +399,13 @@ namespace Opm {
                 GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
                 GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map);
 
+            // cannot be const since it accesses the non-const WellState
+            void gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
+                DeferredLogger& deferred_logger,
+                GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
+                GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
+                GLiftSyncGroups& groups_to_sync);
+
             void extractLegacyCellPvtRegionIndex_();
 
             void extractLegacyDepth_();

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -956,11 +956,9 @@ namespace Opm {
                 for (const auto& well : well_container_) {
                     // NOTE: Only the wells in "group_info" needs to be optimized
                     if (group_info.hasWell(well->name())) {
-                        well->gasLiftOptimizationStage1(
+                        gasLiftOptimizationStage1SingleWell(
-                            this->wellState(), this->groupState(), ebosSimulator_,
+                            well.get(), deferred_logger, prod_wells, glift_wells,
-                            deferred_logger,
+                            group_info, state_map, groups_to_sync
-                            prod_wells, glift_wells, state_map,
-                            group_info, groups_to_sync, this->glift_debug
                         );
                     }
                 }
@@ -1028,6 +1026,35 @@ namespace Opm {
         }
     }
 
+    // NOTE: this method cannot be const since it passes this->wellState()
+    //   (see below) to the GasLiftSingleWell constructor which accepts WellState
+    //   as a non-const reference..
+    template<typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    gasLiftOptimizationStage1SingleWell(WellInterface<TypeTag> *well,
+        DeferredLogger& deferred_logger,
+        GLiftProdWells &prod_wells, GLiftOptWells &glift_wells,
+        GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map,
+        GLiftSyncGroups& sync_groups)
+    {
+        const auto& summary_state = ebosSimulator_.vanguard().summaryState();
+        std::unique_ptr<GasLiftSingleWell> glift
+            = std::make_unique<GasLiftSingleWell>(
+                *well, ebosSimulator_, summary_state,
+                deferred_logger, this->wellState(), this->groupState(),
+                group_info, sync_groups, this->glift_debug);
+        auto state = glift->runOptimize(
+            ebosSimulator_.model().newtonMethod().numIterations());
+        if (state) {
+            state_map.insert({well->name(), std::move(state)});
+            glift_wells.insert({well->name(), std::move(glift)});
+            return;
+        }
+        prod_wells.insert({well->name(), well});
+    }
+
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::

opm/simulators/wells/GasLiftGroupInfo.cpp

@@ -365,6 +365,54 @@ checkNewtonIterationIdxOk_(const std::string &well_name)
     }
 }
 
+void
+GasLiftGroupInfo::
+debugDisplayWellContribution_(
+    const std::string& gr_name, const std::string& well_name,
+    double eff_factor,
+    double well_oil_rate, double well_gas_rate, double well_water_rate,
+    double well_alq,
+    double oil_rate, double gas_rate, double water_rate,
+    double alq
+) const
+{
+    const std::string msg = fmt::format("Group rate for {} : Well {} : "
+        "eff_factor = {}, oil_rate = {}, gas_rate = {}, water_rate = {}, "
+        "alq = {}, New group rates:  oil = {}, gas = {}, water = {}, alq = {}",
+        gr_name, well_name, eff_factor, well_oil_rate, well_gas_rate,
+        well_water_rate, well_alq, oil_rate, gas_rate, water_rate, alq);
+    displayDebugMessage_(msg);
+}
+
+void
+GasLiftGroupInfo::
+debugDisplayUpdatedGroupRates(
+      const std::string& name,
+      double oil_rate, double gas_rate, double water_rate, double alq) const
+{
+
+    const std::string msg = fmt::format("Updated group info for {} : "
+        "oil_rate = {}, gas_rate = {}, water_rate = {}, alq = {}",
+        name, oil_rate, gas_rate, water_rate, alq);
+    displayDebugMessage_(msg);
+}
+
+void
+GasLiftGroupInfo::
+debugEndInitializeGroup(const std::string& name) const
+{
+    const std::string msg = fmt::format("Finished with group {} ...", name);
+    displayDebugMessage_(msg);
+}
+
+void
+GasLiftGroupInfo::
+debugStartInitializeGroup(const std::string& name) const
+{
+    const std::string msg = fmt::format("Initializing group {} ...", name);
+    displayDebugMessage_(msg);
+}
+
 void
 GasLiftGroupInfo::
 displayDebugMessage_(const std::string &msg) const
@@ -417,6 +465,7 @@ GasLiftGroupInfo::
 initializeGroupRatesRecursive_(const Group &group)
 {
     std::array<double,4> rates{};
+    if (this->debug) debugStartInitializeGroup(group.name());
     auto& [oil_rate, water_rate, gas_rate, alq] = rates;
     if (group.wellgroup()) {
         for (const std::string& well_name : group.wells()) {
@@ -439,6 +488,13 @@ initializeGroupRatesRecursive_(const Group &group)
                     gas_rate += (factor * sw_gas_rate);
                     water_rate += (factor * sw_water_rate);
                     alq += (factor * sw_alq);
+                    if (this->debug) {
+                        debugDisplayWellContribution_(
+                            group.name(), well_name, factor,
+                            sw_oil_rate, sw_gas_rate, sw_water_rate, sw_alq,
+                            oil_rate, gas_rate, water_rate, alq
+                        );
+                    }
                 }
             }
         }
@@ -459,6 +515,7 @@ initializeGroupRatesRecursive_(const Group &group)
             alq += (gefac * sg_alq);
         }
     }
+    if (this->debug) debugEndInitializeGroup(group.name());
     std::optional<double> oil_target, gas_target, water_target, liquid_target, max_total_gas, max_alq;
     const auto controls = group.productionControls(this->summary_state_);
     if (group.has_control(Group::ProductionCMode::LRAT)) {
@@ -482,6 +539,10 @@ initializeGroupRatesRecursive_(const Group &group)
         updateGroupIdxMap_(group.name());
         this->group_rate_map_.try_emplace(group.name(),
             oil_rate, gas_rate, water_rate, alq, oil_target, gas_target, water_target, liquid_target, max_total_gas, max_alq);
+        if (this->debug) {
+            debugDisplayUpdatedGroupRates(
+                group.name(), oil_rate, gas_rate, water_rate, alq);
+        }
     }
     return std::make_tuple(oil_rate, gas_rate, water_rate, alq);
 }

opm/simulators/wells/GasLiftGroupInfo.hpp

@@ -114,6 +114,18 @@ public:
 protected:
     bool checkDoGasLiftOptimization_(const std::string& well_name);
     bool checkNewtonIterationIdxOk_(const std::string& well_name);
+    void debugDisplayWellContribution_(
+        const std::string& gr_name, const std::string& well_name,
+        double eff_factor,
+        double well_oil_rate, double well_gas_rate, double well_water_rate,
+        double well_alq,
+        double oil_rate, double gas_rate, double water_rate,
+        double alq
+    ) const;
+    void debugDisplayUpdatedGroupRates(const std::string& name,
+        double oil_rate, double gas_rate, double water_rate, double alq) const;
+    void debugEndInitializeGroup(const std::string& name) const;
+    void debugStartInitializeGroup(const std::string& name) const;
     void displayDebugMessage_(const std::string& msg) const override;
     void displayDebugMessage_(const std::string& msg, const std::string& well_name);
     std::tuple<double, double, double> getProducerWellRates_(const int index);

opm/simulators/wells/GasLiftSingleWell.hpp

@@ -23,16 +23,9 @@
 #include <opm/models/utils/propertysystem.hh>
 #include <opm/models/utils/parametersystem.hh>
 #include <opm/models/discretization/common/fvbaseproperties.hh>
-
-// NOTE: StandardWell.hpp includes ourself (GasLiftSingleWell.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/GasLiftSingleWellGeneric.hpp>
 #include <opm/simulators/wells/GasLiftGroupInfo.hpp>
-
+#include <opm/simulators/wells/WellInterface.hpp>
 
 #include <optional>
 #include <vector>
@@ -45,12 +38,11 @@ namespace Opm
     class GasLiftSingleWell : public GasLiftSingleWellGeneric
     {
         using Simulator = GetPropType<TypeTag, Properties::Simulator>;
-        using StdWell = StandardWell<TypeTag>;
         using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups;
 
     public:
         GasLiftSingleWell(
-            const StdWell &std_well,
+            const WellInterface<TypeTag> &std_well,
             const Simulator &ebos_simulator,
             const SummaryState &summary_state,
             DeferredLogger &deferred_logger,
@@ -69,7 +61,7 @@ namespace Opm
         void setAlqMaxRate_(const GasLiftOpt::Well& well);
 
         const Simulator &ebos_simulator_;
-        const StdWell &std_well_;
+        const WellInterface<TypeTag> &std_well_;
     };
 
 } // namespace Opm

opm/simulators/wells/GasLiftSingleWellGeneric.cpp

@@ -285,7 +285,13 @@ checkThpControl_() const
     const int well_index = this->well_state_.index(this->well_name_).value();
     const Well::ProducerCMode& control_mode =
                          this->well_state_.well(well_index).production_cmode;
-    return control_mode == Well::ProducerCMode::THP;
+    bool thp_control = control_mode == Well::ProducerCMode::THP;
+    if (this->debug) {
+        if (!thp_control) {
+            displayDebugMessage_("Well is not under THP control, skipping iteration..");
+        }
+    }
+    return thp_control;
 }
 
 
@@ -431,6 +437,18 @@ debugShowLimitingTargets_(const LimitedRates& rates) const
     }
 }
 
+void
+GasLiftSingleWellGeneric::
+debugShowProducerControlMode() const
+{
+    const int well_index = this->well_state_.index(this->well_name_).value();
+    const Well::ProducerCMode& control_mode =
+                         this->well_state_.well(well_index).production_cmode;
+    const std::string msg = fmt::format("Current control mode is: {}",
+        Well::ProducerCMode2String(control_mode));
+    displayDebugMessage_(msg);
+}
+
 void
 GasLiftSingleWellGeneric::
 debugShowStartIteration_(double alq, bool increase, double oil_rate)
@@ -1139,6 +1157,7 @@ std::unique_ptr<GasLiftWellState>
 GasLiftSingleWellGeneric::
 runOptimizeLoop_(bool increase)
 {
+    if (this->debug) debugShowProducerControlMode();
     std::unique_ptr<GasLiftWellState> ret_value; // nullptr initially
     auto rates = getInitialRatesWithLimit_();
     if (!rates) return ret_value;
@@ -1160,7 +1179,11 @@ runOptimizeLoop_(bool increase)
     }
 
     OptimizeState state {*this, increase};
-    if (!checkThpControl_()) {
+    auto temp_alq = cur_alq;
+    if (checkThpControl_()) {
+        if (this->debug) debugShowStartIteration_(temp_alq, increase, new_rates.oil);
+    }
+    else {
         // If the well is not under THP control, we can still use the previous
         //   initial adjustment of ALQ by using the well's THP limit to calculate
         //   BHP and then well rates from that.
@@ -1169,8 +1192,6 @@ runOptimizeLoop_(bool increase)
         //   Then gaslift can be reduced while still keeping the group target.
         state.stop_iteration = true;
     }
-    auto temp_alq = cur_alq;
-    if (this->debug) debugShowStartIteration_(temp_alq, increase, new_rates.oil);
     while (!state.stop_iteration && (++state.it <= this->max_iterations_)) {
         if (state.checkRatesViolated(new_rates)) break;
         if (state.checkAlqOutsideLimits(temp_alq, new_rates.oil)) break;
@@ -1524,11 +1545,11 @@ checkAlqOutsideLimits(double alq, [[maybe_unused]] double oil_rate)
     else { // we are decreasing lift gas
         if ( alq == 0 ) {
             ss << "ALQ is zero, cannot decrease further. Stopping iteration.";
-            return true;
+            result = true;
         }
         else if ( alq < 0 ) {
             ss << "Negative ALQ: " << alq << ". Stopping iteration.";
-            return true;
+            result = true;
         }
         // NOTE: A negative min_alq_ means: allocate at least enough lift gas
         //  to enable the well to flow, see WLIFTOPT item 5.

opm/simulators/wells/GasLiftSingleWellGeneric.hpp

@@ -260,6 +260,7 @@ protected:
     void debugShowAlqIncreaseDecreaseCounts_();
     void debugShowBhpAlqTable_();
     void debugShowLimitingTargets_(const LimitedRates& rates) const;
+    void debugShowProducerControlMode() const;
     void debugShowStartIteration_(double alq, bool increase, double oil_rate);
     void debugShowTargets_();
     void displayDebugMessage_(const std::string& msg) const override;

opm/simulators/wells/GasLiftSingleWell_impl.hpp

@@ -21,7 +21,7 @@ namespace Opm {
 
 template<typename TypeTag>
 GasLiftSingleWell<TypeTag>::
-GasLiftSingleWell(const StdWell &std_well,
+GasLiftSingleWell(const WellInterface<TypeTag> &std_well,
                   const Simulator &ebos_simulator,
                   const SummaryState &summary_state,
                   DeferredLogger &deferred_logger,
@@ -167,7 +167,7 @@ setAlqMaxRate_(const GasLiftOpt::Well &well)
         // According to the 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(
+        const auto& table = std_well_.vfpProperties()->getProd()->getTable(
                 this->controls_.vfp_table_number);
         const auto& alq_values = table.getALQAxis();
         // Assume the alq_values are sorted in ascending order, so

opm/simulators/wells/MultisegmentWell.hpp

@@ -52,11 +52,6 @@ namespace Opm
         using typename Base::RateConverterType;
         using typename Base::SparseMatrixAdapter;
         using typename Base::FluidState;
-        using typename Base::GasLiftSingleWell;
-        using typename Base::GLiftProdWells;
-        using typename Base::GLiftOptWells;
-        using typename Base::GLiftWellStateMap;
-        using typename Base::GLiftSyncGroups;
 
         using Base::has_solvent;
         using Base::has_polymer;
@@ -99,21 +94,6 @@ namespace Opm
 
         virtual void initPrimaryVariablesEvaluation() const override;
 
-        virtual void gasLiftOptimizationStage1 (
-            WellState&,
-            const GroupState&,
-            const Simulator&,
-            DeferredLogger&,
-            GLiftProdWells &,
-            GLiftOptWells &,
-            GLiftWellStateMap &,
-            GasLiftGroupInfo &,
-            GLiftSyncGroups &,
-            bool
-        ) const override {
-            // Not implemented yet
-        }
-
         /// updating the well state based the current control mode
         virtual void updateWellStateWithTarget(const Simulator& ebos_simulator,
                                                const GroupState& group_state,
@@ -173,6 +153,12 @@ namespace Opm
                                     double* connII,
                                     DeferredLogger& deferred_logger) const;
 
+        virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
+            const Simulator& ebos_simulator,
+            const SummaryState& summary_state,
+            DeferredLogger& deferred_logger,
+            double alq_value) const override;
+
     protected:
         int number_segments_;
 
@@ -258,19 +244,20 @@ namespace Opm
                                         std::vector<double>& well_flux,
                                         DeferredLogger& deferred_logger) const;
 
-        void computeWellRatesWithBhp(const Simulator& ebosSimulator,
+        virtual void computeWellRatesWithBhp(const Simulator& ebosSimulator,
-                                     const Scalar& bhp,
+                                     const double& bhp,
                                      std::vector<double>& well_flux,
-                                     DeferredLogger& deferred_logger) const;
+                                     DeferredLogger& deferred_logger) const override;
 
         void computeWellRatesWithBhpIterations(const Simulator& ebosSimulator,
                                                const Scalar& bhp,
                                                std::vector<double>& well_flux,
                                                DeferredLogger& deferred_logger) const;
 
-        std::vector<double>
+        std::vector<double> computeWellPotentialWithTHP(
-        computeWellPotentialWithTHP(const Simulator& ebos_simulator,
+                                 const WellState& well_state,
-                                    DeferredLogger& deferred_logger) const;
+                                 const Simulator& ebos_simulator,
+                                 DeferredLogger& deferred_logger) const;
 
         virtual double getRefDensity() const override;
 
@@ -306,9 +293,12 @@ namespace Opm
         bool allDrawDownWrongDirection(const Simulator& ebos_simulator) const;
 
 
-        std::optional<double> computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
+
-                                                       const SummaryState& summary_state,
+        std::optional<double> computeBhpAtThpLimitProd(
-                                                       DeferredLogger& deferred_logger) const;
+            const WellState& well_state,
+            const Simulator& ebos_simulator,
+            const SummaryState& summary_state,
+            DeferredLogger& deferred_logger) const;
 
         std::optional<double> computeBhpAtThpLimitInj(const Simulator& ebos_simulator,
                                                       const SummaryState& summary_state,

opm/simulators/wells/MultisegmentWellGeneric.cpp

@@ -436,11 +436,13 @@ computeBhpAtThpLimitInj(const std::function<std::vector<double>(const double)>&
 template<typename Scalar>
 std::optional<double>
 MultisegmentWellGeneric<Scalar>::
-computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
+computeBhpAtThpLimitProdWithAlq(
-                         const SummaryState& summary_state,
+                   const std::function<std::vector<double>(const double)>& frates,
-                         const double maxPerfPress,
+                   const SummaryState& summary_state,
-                         const double rho,
+                   const double maxPerfPress,
-                         DeferredLogger& deferred_logger) const
+                   const double rho,
+                   DeferredLogger& deferred_logger,
+                   double alq_value) const
 {
     // Given a VFP function returning bhp as a function of phase
     // rates and thp:
@@ -468,10 +470,11 @@ computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>&
     const double vfp_ref_depth = table.getDatumDepth();
     const double thp_limit = baseif_.getTHPConstraint(summary_state);
     const double dp = wellhelpers::computeHydrostaticCorrection(baseif_.refDepth(), vfp_ref_depth, rho, baseif_.gravity());
-    auto fbhp = [this, &controls, thp_limit, dp](const std::vector<double>& rates) {
+    auto fbhp = [this, &controls, thp_limit, dp, alq_value](
+                                  const std::vector<double>& rates) {
         assert(rates.size() == 3);
         return baseif_.vfpProperties()->getProd()
-        ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, controls.alq_value) - dp;
+        ->bhp(controls.vfp_table_number, rates[Water], rates[Oil], rates[Gas], thp_limit, alq_value) - dp;
     };
 
     // Make the flo() function.

opm/simulators/wells/MultisegmentWellGeneric.hpp

@@ -70,11 +70,13 @@ protected:
                                                   const double rho,
                                                   DeferredLogger& deferred_logger) const;
 
-    std::optional<double> computeBhpAtThpLimitProd(const std::function<std::vector<double>(const double)>& frates,
+    std::optional<double> computeBhpAtThpLimitProdWithAlq(
-                                                   const SummaryState& summary_state,
+        const std::function<std::vector<double>(const double)>& frates,
-                                                   const double maxPerfPress,
+        const SummaryState& summary_state,
-                                                   const double rho,
+        const double maxPerfPress,
-                                                   DeferredLogger& deferred_logger) const;
+        const double rho,
+        DeferredLogger& deferred_logger,
+        double alq_value) const;
 
     std::optional<double> bhpMax(const std::function<double(const double)>& fflo,
                                  const double bhp_limit,

opm/simulators/wells/MultisegmentWell_impl.hpp

@@ -300,7 +300,8 @@ namespace Opm
         if (!Base::wellHasTHPConstraints(summaryState) || bhp_controlled_well) {
             computeWellRatesAtBhpLimit(ebosSimulator, well_potentials, deferred_logger);
         } else {
-            well_potentials = computeWellPotentialWithTHP(ebosSimulator, deferred_logger);
+            well_potentials = computeWellPotentialWithTHP(
+                well_state, ebosSimulator, deferred_logger);
         }
         deferred_logger.debug("Cost in iterations of finding well potential for well "
                               + this->name() + ": " + std::to_string(debug_cost_counter_));
@@ -342,7 +343,7 @@ namespace Opm
     void
     MultisegmentWell<TypeTag>::
     computeWellRatesWithBhp(const Simulator& ebosSimulator,
-                            const Scalar& bhp,
+                            const double& bhp,
                             std::vector<double>& well_flux,
                             DeferredLogger& deferred_logger) const
     {
@@ -448,8 +449,10 @@ namespace Opm
     template<typename TypeTag>
     std::vector<double>
     MultisegmentWell<TypeTag>::
-    computeWellPotentialWithTHP(const Simulator& ebos_simulator,
+    computeWellPotentialWithTHP(
-                                DeferredLogger& deferred_logger) const
+          const WellState& well_state,
+          const Simulator& ebos_simulator,
+          DeferredLogger& deferred_logger) const
     {
         std::vector<double> potentials(this->number_of_phases_, 0.0);
         const auto& summary_state = ebos_simulator.vanguard().summaryState();
@@ -472,7 +475,8 @@ namespace Opm
                 computeWellRatesWithBhpIterations(ebos_simulator, bhp, potentials, deferred_logger);
             }
         } else {
-            auto bhp_at_thp_limit = computeBhpAtThpLimitProd(ebos_simulator, summary_state, deferred_logger);
+            auto bhp_at_thp_limit = computeBhpAtThpLimitProd(
+                  well_state, ebos_simulator, summary_state, deferred_logger);
             if (bhp_at_thp_limit) {
                 const auto& controls = well.productionControls(summary_state);
                 const double bhp = std::max(*bhp_at_thp_limit, controls.bhp_limit);
@@ -1293,11 +1297,16 @@ namespace Opm
     template<typename TypeTag>
     void
     MultisegmentWell<TypeTag>::
-    checkOperabilityUnderTHPLimit(const Simulator& ebos_simulator, const WellState& /*well_state*/, DeferredLogger& deferred_logger)
+    checkOperabilityUnderTHPLimit(
+             const Simulator& ebos_simulator,
+             const WellState& well_state,
+             DeferredLogger& deferred_logger)
     {
         const auto& summaryState = ebos_simulator.vanguard().summaryState();
-        const auto obtain_bhp = this->isProducer() ? computeBhpAtThpLimitProd(ebos_simulator, summaryState, deferred_logger)
+        const auto obtain_bhp = this->isProducer()
-        : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
+            ? computeBhpAtThpLimitProd(
+                        well_state, ebos_simulator, summaryState, deferred_logger)
+            : computeBhpAtThpLimitInj(ebos_simulator, summaryState, deferred_logger);
 
         if (obtain_bhp) {
             this->operability_status_.can_obtain_bhp_with_thp_limit = true;
@@ -1734,15 +1743,30 @@ namespace Opm
     }
 
 
+    template<typename TypeTag>
+    std::optional<double>
+    MultisegmentWell<TypeTag>::
+    computeBhpAtThpLimitProd(const WellState& well_state,
+                             const Simulator& ebos_simulator,
+                             const SummaryState& summary_state,
+                             DeferredLogger& deferred_logger) const
+    {
+        return this->MultisegmentWell<TypeTag>::computeBhpAtThpLimitProdWithAlq(
+                                               ebos_simulator,
+                                               summary_state,
+                                               deferred_logger,
+                                               this->getALQ(well_state));
+    }
 
 
 
     template<typename TypeTag>
     std::optional<double>
     MultisegmentWell<TypeTag>::
-    computeBhpAtThpLimitProd(const Simulator& ebos_simulator,
+    computeBhpAtThpLimitProdWithAlq(const Simulator& ebos_simulator,
-                             const SummaryState& summary_state,
+                                    const SummaryState& summary_state,
-                             DeferredLogger& deferred_logger) const
+                                    DeferredLogger& deferred_logger,
+                                    double alq_value) const
     {
         // Make the frates() function.
         auto frates = [this, &ebos_simulator, &deferred_logger](const double bhp) {
@@ -1757,11 +1781,12 @@ namespace Opm
         };
 
         auto bhpAtLimit = this->MultisegmentWellGeneric<Scalar>::
-               computeBhpAtThpLimitProd(frates,
+               computeBhpAtThpLimitProdWithAlq(frates,
                                         summary_state,
                                         maxPerfPress(ebos_simulator),
                                         getRefDensity(),
-                                        deferred_logger);
+                                        deferred_logger,
+                                        alq_value);
 
        if(bhpAtLimit)
            return bhpAtLimit;
@@ -1776,11 +1801,12 @@ namespace Opm
        };
 
        return this->MultisegmentWellGeneric<Scalar>::
-              computeBhpAtThpLimitProd(fratesIter,
+              computeBhpAtThpLimitProdWithAlq(fratesIter,
                                        summary_state,
                                        maxPerfPress(ebos_simulator),
                                        getRefDensity(),
-                                       deferred_logger);
+                                       deferred_logger,
+                                       alq_value);
 
     }
 

opm/simulators/wells/StandardWell.hpp

@@ -31,8 +31,6 @@
 #include <opm/simulators/wells/WellInterface.hpp>
 #include <opm/simulators/wells/WellProdIndexCalculator.hpp>
 #include <opm/simulators/wells/ParallelWellInfo.hpp>
-#include <opm/simulators/wells/GasLiftSingleWell.hpp>
-#include <opm/simulators/wells/GasLiftGroupInfo.hpp>
 
 #include <opm/models/blackoil/blackoilpolymermodules.hh>
 #include <opm/models/blackoil/blackoilsolventmodules.hh>
@@ -83,11 +81,6 @@ namespace Opm
         using typename Base::SparseMatrixAdapter;
         using typename Base::FluidState;
         using typename Base::RateVector;
-        using typename Base::GasLiftSingleWell;
-        using typename Base::GLiftOptWells;
-        using typename Base::GLiftProdWells;
-        using typename Base::GLiftWellStateMap;
-        using typename Base::GLiftSyncGroups;
 
         using Base::has_solvent;
         using Base::has_zFraction;
@@ -202,19 +195,6 @@ namespace Opm
             return this->param_.matrix_add_well_contributions_;
         }
 
-        virtual void gasLiftOptimizationStage1 (
-            WellState& well_state,
-            const GroupState& group_state,
-            const Simulator& ebosSimulator,
-            DeferredLogger& deferred_logger,
-            GLiftProdWells &prod_wells,
-            GLiftOptWells &glift_wells,
-            GLiftWellStateMap &state_map,
-            GasLiftGroupInfo &group_info,
-            GLiftSyncGroups &sync_groups,
-            bool glift_debug
-        ) const override;
-
         /* returns BHP */
         double computeWellRatesAndBhpWithThpAlqProd(const Simulator &ebos_simulator,
                                const SummaryState &summary_state,
@@ -229,19 +209,17 @@ namespace Opm
             std::vector<double> &potentials,
             double alq) const;
 
-        // NOTE: Cannot be protected since it is used by GasLiftRuntime
+        virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
-        std::optional<double> computeBhpAtThpLimitProdWithAlq(
             const Simulator& ebos_simulator,
             const SummaryState& summary_state,
             DeferredLogger& deferred_logger,
-            double alq_value) const;
+            double alq_value) const override;
 
-        // NOTE: Cannot be protected since it is used by GasLiftRuntime
+        virtual void computeWellRatesWithBhp(
-        void computeWellRatesWithBhp(
             const Simulator& ebosSimulator,
             const double& bhp,
             std::vector<double>& well_flux,
-            DeferredLogger& deferred_logger) const;
+            DeferredLogger& deferred_logger) const override;
 
         // NOTE: These cannot be protected since they are used by GasLiftRuntime
         using Base::phaseUsage;

opm/simulators/wells/StandardWell_impl.hpp

@@ -1812,38 +1812,6 @@ namespace Opm
                                              alq);
     }
 
-    template<typename TypeTag>
-    void
-    StandardWell<TypeTag>::
-    gasLiftOptimizationStage1(
-                       WellState& well_state,
-                       const GroupState& group_state,
-                       const Simulator& ebos_simulator,
-                       DeferredLogger& deferred_logger,
-                       GLiftProdWells &prod_wells,
-                       GLiftOptWells &glift_wells,
-                       GLiftWellStateMap &glift_state_map,
-                       GasLiftGroupInfo &group_info,
-                       GLiftSyncGroups &sync_groups,
-                       bool glift_debug
-    ) const
-    {
-        const auto& summary_state = ebos_simulator.vanguard().summaryState();
-        std::unique_ptr<GasLiftSingleWell> glift
-            = std::make_unique<GasLiftSingleWell>(
-                *this, ebos_simulator, summary_state,
-                deferred_logger, well_state, group_state, group_info,
-                sync_groups, glift_debug);
-        auto state = glift->runOptimize(
-            ebos_simulator.model().newtonMethod().numIterations());
-        if (state) {
-            glift_state_map.insert({this->name(), std::move(state)});
-            glift_wells.insert({this->name(), std::move(glift)});
-            return;
-        }
-        prod_wells.insert({this->name(), this});
-    }
-
     template<typename TypeTag>
     void
     StandardWell<TypeTag>::

opm/simulators/wells/WellInterface.hpp

@@ -57,6 +57,7 @@ namespace Opm {
 #include <opm/material/densead/Evaluation.hpp>
 
 #include <opm/simulators/wells/WellInterfaceIndices.hpp>
+#include <opm/simulators/timestepping/ConvergenceReport.hpp>
 
 #include <cassert>
 #include <vector>
@@ -159,17 +160,18 @@ public:
                         const GroupState& group_state,
                         DeferredLogger& deferred_logger);
 
-    virtual void gasLiftOptimizationStage1 (
+    virtual void computeWellRatesWithBhp(
-        WellState& well_state,
-        const GroupState& group_state,
         const Simulator& ebosSimulator,
+        const double& bhp,
+        std::vector<double>& well_flux,
+        DeferredLogger& deferred_logger
+    ) const = 0;
+
+    virtual std::optional<double> computeBhpAtThpLimitProdWithAlq(
+        const Simulator& ebos_simulator,
+        const SummaryState& summary_state,
         DeferredLogger& deferred_logger,
-        GLiftProdWells& prod_wells,
+        double alq_value
-        GLiftOptWells& glift_wells,
-        GLiftWellStateMap& state_map,
-        GasLiftGroupInfo &group_info,
-        GLiftSyncGroups &sync_groups,
-        bool glift_debug
     ) const = 0;
 
     /// using the solution x to recover the solution xw for wells and applying