Cleanup code in GasLiftSingleWellGeneric

Introduces two new data types BasicRates and LimitedRates to capture
oil, gas, and water rates, and whether they have been limited by well
or group targets. This reduces the number of variables that are passed
to and returned from various methods and thus makes the code easier to
read.

opm/simulators/wells/GasLiftGroupInfo.cpp

@@ -102,7 +102,6 @@ getRate(Rate rate_type, const std::string& group_name) const
     }
 }
 
-
 std::tuple<double, double, double, double>
 GasLiftGroupInfo::
 getRates(const int group_idx) const
@@ -160,6 +159,14 @@ groupIdxToName(int group_idx) const
     return *group_name;
 }
 
+bool
+GasLiftGroupInfo::
+hasAnyTarget(const std::string& group_name) const
+{
+    return oilTarget(group_name)   || gasTarget(group_name)
+        || waterTarget(group_name) || liquidTarget(group_name);
+}
+
 bool
 GasLiftGroupInfo::
 hasWell(const std::string& well_name)

opm/simulators/wells/GasLiftGroupInfo.hpp

@@ -96,6 +96,7 @@ public:
     std::optional<double> getTarget(
         Rate rate_type, const std::string& group_name) const;
     const std::string& groupIdxToName(int group_idx) const;
+    bool hasAnyTarget(const std::string& group_name) const;
     bool hasWell(const std::string& well_name);
     void initialize();
     std::optional<double> liquidTarget(const std::string& group_name) const;
@@ -110,9 +111,9 @@ public:
     void updateRate(int idx, double oil_rate, double gas_rate, double water_rate, double alq);
     const Well2GroupMap& wellGroupMap() { return well_group_map_; }
 protected:
-    void displayDebugMessage_(const std::string& msg) const override;
     bool checkDoGasLiftOptimization_(const std::string& well_name);
     bool checkNewtonIterationIdxOk_(const std::string& well_name);
+    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);
     std::tuple<double, double, double, double>

opm/simulators/wells/GasLiftSingleWell.hpp

@@ -64,9 +64,8 @@ namespace Opm
 
     private:
         std::optional<double> computeBhpAtThpLimit_(double alq) const override;
-        void computeWellRates_(
+        BasicRates computeWellRates_(
-            double bhp, std::vector<double> &potentials, bool debug_output=true) const override;
+            double bhp, bool bhp_is_limited, bool debug_output=true) const override;
-
         void setAlqMaxRate_(const GasLiftOpt::Well& well);
 
         const Simulator &ebos_simulator_;

opm/simulators/wells/GasLiftSingleWellGeneric.hpp

@@ -107,6 +107,101 @@ protected:
         bool glift_debug
     );
 
+    struct LimitedRates;
+    struct BasicRates
+    {
+        BasicRates(const BasicRates& rates) :
+            oil{rates.oil},
+            gas{rates.gas},
+            water{rates.water},
+            bhp_is_limited{rates.bhp_is_limited}
+        {}
+        BasicRates(double oil_, double gas_, double water_, bool bhp_is_limited_) :
+            oil{oil_},
+            gas{gas_},
+            water{water_},
+            bhp_is_limited{bhp_is_limited_}
+        {}
+        BasicRates& operator=(const BasicRates& rates) {
+            oil = rates.oil;
+            gas = rates.gas;
+            water = rates.water;
+            bhp_is_limited = rates.bhp_is_limited;
+            return *this;
+        }
+        // This copy constructor cannot be defined inline here since LimitedRates
+        //   has not been defined yet (it is defined below). Instead it is defined in
+        //   in the .cpp file
+        BasicRates(const LimitedRates& rates);
+        double operator[](Rate rate_type) const {
+            switch (rate_type) {
+            case Rate::oil:
+                return this->oil;
+            case Rate::gas:
+                return this->gas;
+            case Rate::water:
+                return this->water;
+            case Rate::liquid:
+                return this->oil + this->water;
+            default:
+                throw std::runtime_error("This should not happen");
+            }
+        }
+
+        double oil, gas, water;
+        bool bhp_is_limited;
+    };
+
+    struct LimitedRates : public BasicRates
+    {
+        enum class LimitType {well, group, none};
+        LimitedRates(
+              double oil_, double gas_, double water_,
+              bool oil_is_limited_, bool gas_is_limited_,
+              bool water_is_limited_, bool bhp_is_limited_,
+              std::optional<Rate> oil_limiting_target_,
+              std::optional<Rate> water_limiting_target_
+        ) :
+            BasicRates(oil_, gas_, water_, bhp_is_limited_),
+            oil_is_limited{oil_is_limited_},
+            gas_is_limited{gas_is_limited_},
+            water_is_limited{water_is_limited_},
+            oil_limiting_target{oil_limiting_target_},
+            water_limiting_target{water_limiting_target_}
+        {
+            set_initial_limit_type_();
+        }
+
+        LimitedRates(
+              const BasicRates& rates,
+              bool oil_is_limited_, bool gas_is_limited_,
+              bool water_is_limited_
+        ) :
+            BasicRates(rates),
+            oil_is_limited{oil_is_limited_},
+            gas_is_limited{gas_is_limited_},
+            water_is_limited{water_is_limited_}
+        {
+            set_initial_limit_type_();
+        }
+
+        bool limited() const {
+            return oil_is_limited || gas_is_limited || water_is_limited;
+        }
+        // For a given ALQ value, were the rates limited due to group targets
+        //   or due to well targets?
+        LimitType limit_type;
+        bool oil_is_limited;
+        bool gas_is_limited;
+        bool water_is_limited;
+        std::optional<Rate> oil_limiting_target;
+        std::optional<Rate> water_limiting_target;
+    private:
+        void set_initial_limit_type_() {
+            limit_type = limited() ? LimitType::well : LimitType::none;
+        }
+    };
+
     struct OptimizeState
     {
         OptimizeState( GasLiftSingleWellGeneric& parent_, bool increase_ ) :
@@ -129,126 +224,91 @@ protected:
         bool checkAlqOutsideLimits(double alq, double oil_rate);
         bool checkEcoGradient(double gradient);
         bool checkGroupALQrateExceeded(double delta_alq);
-        bool checkNegativeOilRate(double oil_rate);
-        bool checkThpControl();
         bool checkOilRateExceedsTarget(double oil_rate);
-        bool checkRate(double rate, double limit, const std::string &rate_str) const;
+        bool checkRatesViolated(const LimitedRates& rates) const;
-        bool checkWellRatesViolated(std::vector<double> &potentials);
-        bool computeBhpAtThpLimit(double alq);
         void debugShowIterationInfo(double alq);
         double getBhpWithLimit();
-        void updateGroupRates(double delta_oil, double delta_gas, double delta_water, double delta_alq);
         void warn_(std::string msg) {parent.displayWarning_(msg);}
     };
 
-    std::pair<std::optional<double>, bool>
+    std::pair<std::optional<double>, bool> addOrSubtractAlqIncrement_(
-    addOrSubtractAlqIncrement_(double alq, bool increase) const;
+                            double alq, bool increase) const;
-
     double calcEcoGradient_(double oil_rate, double new_oil_rate,
                             double gas_rate, double new_gas_rate, bool increase) const;
-
     bool checkALQequal_(double alq1, double alq2) const;
     bool checkGroupTargetsViolated(
-        double delta_oil, double delta_gas, double delta_water) const;
+                      const BasicRates& rates, const BasicRates& new_rates) const;
     bool checkInitialALQmodified_(double alq, double initial_alq) const;
-
+    bool checkThpControl_() const;
-    bool checkWellRatesViolated_(std::vector<double>& potentials,
-                                 const std::function<bool(double, double, const std::string &)>& callback,
-                                 bool increase);
-
     virtual std::optional<double> computeBhpAtThpLimit_(double alq) const = 0;
-    virtual void computeWellRates_(double bhp,
+    std::optional<BasicRates> computeInitialWellRates_() const;
-                                   std::vector<double>& potentials,
+    std::optional<LimitedRates> computeLimitedWellRatesWithALQ_(double alq) const;
-                                   bool debug_output = true) const = 0;
+    virtual BasicRates computeWellRates_(double bhp, bool bhp_is_limited,                                                             bool debug_output = true) const = 0;
-
+    std::optional<BasicRates> computeWellRatesWithALQ_(double alq) const;
-    bool computeInitialWellRates_(std::vector<double>& potentials);
-
     void debugCheckNegativeGradient_(double grad, double alq, double new_alq,
                                      double oil_rate, double new_oil_rate,
                                      double gas_rate, double new_gas_rate,
                                      bool increase) const;
-
     void debugShowAlqIncreaseDecreaseCounts_();
-
     void debugShowBhpAlqTable_();
-
+    void debugShowLimitingTargets_(const LimitedRates& rates) const;
     void debugShowStartIteration_(double alq, bool increase, double oil_rate);
-
     void debugShowTargets_();
-
     void displayDebugMessage_(const std::string& msg) const override;
     void displayWarning_(const std::string& warning);
-
     std::pair<double, bool> getBhpWithLimit_(double bhp) const;
     std::pair<double, bool> getGasRateWithLimit_(
-                           const std::vector<double>& potentials) const;
+                           const BasicRates& rates) const;
     std::pair<double, bool> getGasRateWithGroupLimit_(
                            double new_gas_rate, double gas_rate) const;
-    std::tuple<double,double,double,bool,bool,bool> getInitialRatesWithLimit_(
+    std::optional<LimitedRates> getInitialRatesWithLimit_() const;
-                           const std::vector<double>& potentials);
+    LimitedRates getLimitedRatesFromRates_(const BasicRates& rates) const;
     std::tuple<double,double,bool,bool> getLiquidRateWithGroupLimit_(
                            const double new_oil_rate, const double oil_rate,
                            const double new_water_rate, const double water_rate) const;
     std::pair<double, bool> getOilRateWithGroupLimit_(
                            double new_oil_rate, double oil_rate) const;
-    std::pair<double, bool> getOilRateWithLimit_(
+    std::pair<double, bool> getOilRateWithLimit_(const BasicRates& rates) const;
-                           const std::vector<double>& potentials) const;
+    std::pair<double, std::optional<Rate>> getOilRateWithLimit2_(
+                           const BasicRates& rates) const;
     double getProductionTarget_(Rate rate) const;
-    double getRate_(Rate rate_type, const std::vector<double>& potentials) const;
+    double getRate_(Rate rate_type, const BasicRates& rates) const;
-    std::pair<double, bool> getRateWithLimit_(
+    std::pair<double, std::optional<Rate>> getRateWithLimit_(
-                           Rate rate_type, const std::vector<double>& potentials) const;
+                    Rate rate_type, const BasicRates& rates) const;
     std::tuple<double, const std::string*, double> getRateWithGroupLimit_(
                  Rate rate_type, const double new_rate, const double old_rate) const;
     std::pair<double, bool> getWaterRateWithGroupLimit_(
                            double new_water_rate, double water_rate) const;
-    std::pair<double, bool> getWaterRateWithLimit_(
+    std::pair<double, bool> getWaterRateWithLimit_(const BasicRates& rates) const;
-                           const std::vector<double>& potentials) const;
+    std::pair<double, std::optional<Rate>> getWaterRateWithLimit2_(
+                           const BasicRates& rates) const;
     bool hasProductionControl_(Rate rate) const;
-
+    std::pair<LimitedRates, double> increaseALQtoPositiveOilRate_(
-    std::tuple<double,double,bool,bool,double>
+                           double alq, const LimitedRates& orig_rates) const;
-    increaseALQtoPositiveOilRate_(double alq,
+    std::pair<LimitedRates, double> increaseALQtoMinALQ_(
-                                  double oil_rate,
+                           double alq, const LimitedRates& orig_rates) const;
-                                  double gas_rate,
-                                  bool oil_is_limited,
-                                  bool gas_is_limited,
-                                  std::vector<double>& potentials);
-
-    std::tuple<double,double,bool,bool,double>
-    increaseALQtoMinALQ_(double alq,
-                         double oil_rate,
-                         double gas_rate,
-                         bool oil_is_limited,
-                         bool gas_is_limited,
-                         std::vector<double>& potentials);
-
     void logSuccess_(double alq,
                      const int iteration_idx);
-    std::tuple<double,double,double, double,bool,bool,bool>
+    std::pair<LimitedRates, double> maybeAdjustALQbeforeOptimizeLoop_(
-      maybeAdjustALQbeforeOptimizeLoop_(
+                        const LimitedRates& rates, double alq, bool increase) const;
-          bool increase, double alq, double oil_rate, double gas_rate, double water_rate,
+    std::pair<LimitedRates, double> reduceALQtoGroupTarget(
-          bool oil_is_limited, bool gas_is_limited, bool water_is_limited, std::vector<double> &potentials);
+                        double alq, const LimitedRates& rates) const;
-    std::tuple<double,double,double,double>
+    std::pair<LimitedRates, double> reduceALQtoWellTarget_(
-      reduceALQtoGroupTarget(double alq, double oil_rate, double gas_rate,
+                        double alq, const LimitedRates& rates) const;
-                             double water_rate, std::vector<double> &potentials) const;
-    std::tuple<double,double,double, bool, bool,bool,double>
-      reduceALQtoWellTarget_(double alq, double oil_rate, double gas_rate,
-                             double water_rate, bool oil_is_limited,
-                             bool gas_is_limited, bool water_is_limited,
-                             std::vector<double> &potentials);
-
     std::unique_ptr<GasLiftWellState> runOptimize1_();
     std::unique_ptr<GasLiftWellState> runOptimize2_();
     std::unique_ptr<GasLiftWellState> runOptimizeLoop_(bool increase);
-
     void setAlqMinRate_(const GasLiftOpt::Well& well);
-
     std::unique_ptr<GasLiftWellState> tryIncreaseLiftGas_();
     std::unique_ptr<GasLiftWellState> tryDecreaseLiftGas_();
-
+    void updateGroupRates_(
+        const LimitedRates& rates,
+        const LimitedRates& new_rates,
+        double delta_alq) const;
+    LimitedRates updateRatesToGroupLimits_(
+        const LimitedRates& rates, const LimitedRates& new_rates) const;
     void updateWellStateAlqFixedValue_(const GasLiftOpt::Well& well);
-
     bool useFixedAlq_(const GasLiftOpt::Well& well);
-
     void warnMaxIterationsExceeded_();
 
     const GroupState& group_state_;

opm/simulators/wells/GasLiftSingleWell_impl.hpp

@@ -97,14 +97,11 @@ GasLiftSingleWell(const StdWell &std_well,
  ****************************************/
 
 template<typename TypeTag>
-void
+GasLiftSingleWellGeneric::BasicRates
 GasLiftSingleWell<TypeTag>::
-computeWellRates_(
+computeWellRates_( double bhp, bool bhp_is_limited, bool debug_output ) const
-    double bhp, std::vector<double> &potentials, bool debug_output) const
 {
-    // NOTE: If we do not clear the potentials here, it will accumulate
+    std::vector<double> potentials(this->num_phases_, 0.0);
-    //   the new potentials to the old values..
-    std::fill(potentials.begin(), potentials.end(), 0.0);
     this->std_well_.computeWellRatesWithBhp(
         this->ebos_simulator_, bhp, potentials, this->deferred_logger_);
     if (debug_output) {
@@ -118,6 +115,11 @@ computeWellRates_(
     for (auto& potential : potentials) {
         potential = std::min(0.0, potential);
     }
+    return {-potentials[this->oil_pos_],
+            -potentials[this->gas_pos_],
+            -potentials[this->water_pos_],
+            bhp_is_limited
+    };
 }
 
 template<typename TypeTag>

opm/simulators/wells/GasLiftStage2.cpp

@@ -128,7 +128,11 @@ calcIncOrDecGrad_(
     // only applies to wells in the well_state_map (i.e. wells on this rank)
     if(this->well_state_map_.count(well_name) == 0)
         return std::nullopt;
-
+    if (this->debug) {
+        const std::string msg = fmt::format("well {} : calculating {} gradient..",
+            well_name, (increase ? "incremental" : "decremental"));
+        displayDebugMessage_(msg);
+    }
     GasLiftWellState &state = *(this->well_state_map_.at(well_name).get());
     if (checkRateAlreadyLimited_(state, increase)) {
         /*
@@ -191,7 +195,6 @@ checkRateAlreadyLimited_(GasLiftWellState &state, bool increase)
     return false;
 }
 
-
 GasLiftStage2::GradInfo
 GasLiftStage2::
 deleteGrad_(const std::string &name, bool increase)