Merge pull request #4573 from hakonhagland/glift_grad

Fix bug when recalculating gaslift gradient

opm/simulators/wells/GasLiftSingleWell.hpp

@@ -55,7 +55,7 @@ namespace Opm
         const WellInterfaceGeneric &getWell() const override { return well_; }
 
     private:
-        std::optional<double> computeBhpAtThpLimit_(double alq) const override;
+        std::optional<double> computeBhpAtThpLimit_(double alq, bool debug_ouput=true) const override;
         BasicRates computeWellRates_(
             double bhp, bool bhp_is_limited, bool debug_output=true) const override;
         void setAlqMaxRate_(const GasLiftWell& well);

opm/simulators/wells/GasLiftSingleWellGeneric.cpp

@@ -85,8 +85,14 @@ GasLiftSingleWellGeneric::GasLiftSingleWellGeneric(DeferredLogger& deferred_logg
  ****************************************/
 // NOTE: Used from GasLiftStage2
 std::optional<GasLiftSingleWellGeneric::GradInfo>
-GasLiftSingleWellGeneric::calcIncOrDecGradient(
+GasLiftSingleWellGeneric::calcIncOrDecGradient( double oil_rate,
-    double oil_rate, double gas_rate, double water_rate, double alq, const std::string& gr_name_dont_limit, bool increase) const
+                                                double gas_rate,
+                                                double water_rate,
+                                                double alq,
+                                                const std::string& gr_name_dont_limit,
+                                                bool increase,
+                                                bool debug_output
+    ) const
 {
     auto [new_alq_opt, alq_is_limited] = addOrSubtractAlqIncrement_(alq, increase);
     // TODO: What to do if ALQ is limited and new_alq != alq?
@@ -99,15 +105,10 @@ GasLiftSingleWellGeneric::calcIncOrDecGradient(
     if (checkGroupALQrateExceeded(delta_alq, gr_name_dont_limit))
         return std::nullopt;
 
-    if (auto bhp = computeBhpAtThpLimit_(new_alq)) {
+    if (auto bhp = computeBhpAtThpLimit_(new_alq, debug_output)) {
         auto [new_bhp, bhp_is_limited] = getBhpWithLimit_(*bhp);
         // TODO: What to do if BHP is limited?
-        auto rates = computeWellRates_(new_bhp, bhp_is_limited);
+        auto rates = computeWellRates_(new_bhp, bhp_is_limited, debug_output);
-        // double new_oil_rate, new_gas_rate, new_water_rate;
-        // bool oil_is_limited, gas_is_limited, water_is_limited;
-        // std::tie(new_oil_rate, oil_is_limited) = getOilRateWithLimit_(rates);
-        // std::tie(new_gas_rate, gas_is_limited) = getGasRateWithLimit_(rates);
-        // std::tie(new_water_rate, water_is_limited) = getWaterRateWithLimit_(rates);
         const auto ratesLimited = getLimitedRatesFromRates_(rates);
         BasicRates oldrates = {oil_rate, gas_rate, water_rate, false};
         const auto new_rates = updateRatesToGroupLimits_(oldrates, ratesLimited, gr_name_dont_limit);
@@ -288,7 +289,7 @@ computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const
     auto alq = this->orig_alq_;
     double new_alq = alq;
     std::optional<double> bhp;
-    while (alq <= (this->max_alq_ + this->increment_)) {
+    while ((alq < this->max_alq_) || checkALQequal_(alq, this->max_alq_)) {
         if (bhp = computeBhpAtThpLimit_(alq); bhp) {
             new_alq = alq;
             break;
@@ -321,7 +322,8 @@ GasLiftSingleWellGeneric::computeInitialWellRates_() const
                                                 rates->water);
             displayDebugMessage_(msg);
         }
-    } else {
+    }
+    else {
         displayDebugMessage_("Aborting optimization.");
     }
     return {rates, initial_alq};
@@ -1531,10 +1533,10 @@ GasLiftSingleWellGeneric::OptimizeState::checkAlqOutsideLimits(double alq, [[may
         }
     } else { // we are decreasing lift gas
         if (alq == 0) {
-            ss << "ALQ is zero, cannot decrease further. Stopping iteration.";
+            ss << "ALQ is zero, cannot decrease further. Stopping iteration. ";
             result = true;
         } else if (alq < 0) {
-            ss << "Negative ALQ: " << alq << ". Stopping iteration.";
+            ss << "Negative ALQ: " << alq << ". Stopping iteration. ";
             result = true;
         }
         // NOTE: A negative min_alq_ means: allocate at least enough lift gas
@@ -1593,7 +1595,7 @@ GasLiftSingleWellGeneric::checkGroupALQrateExceeded(double delta_alq, const std:
     const auto& pairs = group_info_.getWellGroups(well_name_);
     for (const auto& [group_name, efficiency] : pairs) {
         // in stage 2 we don't want to limit the rate to the group
-        // target we are trying to redistribute the gaslift within 
+        // target we are trying to redistribute the gaslift within
         if (gr_name_dont_limit == group_name)
             continue;
         auto max_alq_opt = group_info_.maxAlq(group_name);

opm/simulators/wells/GasLiftSingleWellGeneric.hpp

@@ -95,9 +95,11 @@ public:
 
     std::optional<GradInfo> calcIncOrDecGradient(double oil_rate, double gas_rate,
                                                  double water_rate,
-                                                 double alq, 
+                                                 double alq,
                                                  const std::string& gr_name_dont_limit,
-                                                 bool increase) const;
+                                                 bool increase,
+                                                 bool debug_output = true
+                                                 ) const;
 
     std::unique_ptr<GasLiftWellState> runOptimize(const int iteration_idx);
 
@@ -253,11 +255,11 @@ protected:
                       const BasicRates& rates, const BasicRates& new_rates) const;
     bool checkInitialALQmodified_(double alq, double initial_alq) const;
     virtual bool checkThpControl_() const = 0;
-    virtual std::optional<double> computeBhpAtThpLimit_(double alq) const = 0;
+    virtual std::optional<double> computeBhpAtThpLimit_(double alq, bool debug_output = true) const = 0;
     std::pair<std::optional<double>,double> computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const;
     std::pair<std::optional<BasicRates>,double> computeInitialWellRates_() const;
     std::optional<LimitedRates> computeLimitedWellRatesWithALQ_(double alq) const;
-    virtual BasicRates computeWellRates_(double bhp, bool bhp_is_limited,                                                             bool debug_output = true) const = 0;
+    virtual BasicRates computeWellRates_(double bhp, bool bhp_is_limited, bool debug_output = true) const = 0;
     std::optional<BasicRates> computeWellRatesWithALQ_(double alq) const;
     void debugCheckNegativeGradient_(double grad, double alq, double new_alq,
                                      double oil_rate, double new_oil_rate,

opm/simulators/wells/GasLiftSingleWell_impl.hpp

@@ -128,7 +128,7 @@ computeWellRates_( double bhp, bool bhp_is_limited, bool debug_output ) const
 template<typename TypeTag>
 std::optional<double>
 GasLiftSingleWell<TypeTag>::
-computeBhpAtThpLimit_(double alq) const
+computeBhpAtThpLimit_(double alq, bool debug_output) const
 {
     auto bhp_at_thp_limit = this->well_.computeBhpAtThpLimitProdWithAlq(
         this->ebos_simulator_,
@@ -137,11 +137,14 @@ computeBhpAtThpLimit_(double alq) const
         this->deferred_logger_);
     if (bhp_at_thp_limit) {
         if (*bhp_at_thp_limit < this->controls_.bhp_limit) {
-            const std::string msg = fmt::format(
+            if (debug_output && this->debug) {
-                "Computed bhp ({}) from thp limit is below bhp limit ({}), (ALQ = {})."
+                const std::string msg = fmt::format(
-                " Using bhp limit instead",
+                    "Computed bhp ({}) from thp limit is below bhp limit ({}), (ALQ = {})."
-                *bhp_at_thp_limit, this->controls_.bhp_limit, alq);
+                    " Using bhp limit instead",
-            displayDebugMessage_(msg);
+                    *bhp_at_thp_limit, this->controls_.bhp_limit, alq
+                );
+                displayDebugMessage_(msg);
+            }
             bhp_at_thp_limit = this->controls_.bhp_limit;
         }
         //bhp_at_thp_limit = std::max(*bhp_at_thp_limit, this->controls_.bhp_limit);

opm/simulators/wells/GasLiftStage2.cpp

@@ -96,6 +96,8 @@ runOptimize()
 // Update GasLiftWellState and WellState for "well_name" to the
 //   new ALQ value and related data (the data has already been computed and
 //   saved in "grad_map")
+// INPUT: grad_map : map of incremental (if "add" is true) or decremental
+//        (if "add" is false) GradInfo structs for each well name.
 void
 GasLiftStage2::
 addOrRemoveALQincrement_(GradMap &grad_map, const std::string& well_name, bool add)
@@ -141,35 +143,35 @@ 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)) {
+    if (checkRateAlreadyLimited_(well_name, state, increase)) {
-        /*
-        const std::string msg = fmt::format(
-            "well {} : not able to obtain {} gradient",
-            well_name,
-            (increase ? "incremental" : "decremental")
-        );
-        displayDebugMessage_(msg);
-        */
         return std::nullopt;
     }
     else {
         auto [oil_rate, gas_rate] = state.getRates();
         auto alq = state.alq();
-        auto grad = gs_well.calcIncOrDecGradient(oil_rate, gas_rate, state.waterRate(), alq, gr_name_dont_limit, increase);
+        auto grad = gs_well.calcIncOrDecGradient(
-        if (grad) {
+            oil_rate, gas_rate, state.waterRate(), alq, gr_name_dont_limit, increase, /*debug_output=*/false);
-            const std::string msg = fmt::format(
+        if (this->debug) {
-              "well {} : adding {} gradient = {}",
+            if (grad) {
-              well_name,
+                const std::string msg = fmt::format(
-              (increase ? "incremental" : "decremental"),
+                    "well {} : alq = {} : adding {} gradient = {}",
-              grad->grad
+                    well_name,
-            );
+                    alq,
-            displayDebugMessage_(msg);
+                    (increase ? "incremental" : "decremental"),
+                    grad->grad
+                );
+                displayDebugMessage_(msg);
+            }
+            else {
+                const std::string msg = fmt::format(
+                    "well {} : alq = {} : failed to compute {} gradient",
+                    well_name,
+                    alq,
+                    (increase ? "incremental" : "decremental")
+                );
+                displayDebugMessage_(msg);
+            }
         }
         return grad;
     }
@@ -177,7 +179,7 @@ calcIncOrDecGrad_(
 
 bool
 GasLiftStage2::
-checkRateAlreadyLimited_(GasLiftWellState &state, bool increase)
+checkRateAlreadyLimited_(const std::string& well_name, GasLiftWellState &state, bool increase)
 {
     auto current_increase = state.increase();
     bool do_check = false;
@@ -197,10 +199,12 @@ checkRateAlreadyLimited_(GasLiftWellState &state, bool increase)
     if (do_check) {
         if (state.gasIsLimited() || state.oilIsLimited() || state.alqIsLimited() || state.waterIsLimited()) {
             const std::string msg = fmt::format(
-                "{} gradient : skipping since {} was limited in previous step",
+                "Well {} : alq = {} : skipping {} gradient since {} was limited in previous step",
+                well_name,
+                state.alq(),
                 (increase ? "incremental" : "decremental"),
-                (state.oilIsLimited() ? "oil" :
+                (state.oilIsLimited() ? "oil" : (state.gasIsLimited() ? "gas" : "alq"))
-                    (state.gasIsLimited() ? "gas" : "alq")));
+            );
             displayDebugMessage_(msg);
             return true;
         }
@@ -405,14 +409,12 @@ optimizeGroup_(const Group &group)
 {
     const auto& group_name = group.name();
     const auto prod_control = this->group_state_.production_control(group_name);
-    //if (group.has_control(Group::ProductionCMode::ORAT) || group.has_control(Group::ProductionCMode::LRAT)
-    //                   || max_glift || max_total_gas)
-    // NOTE: it only makes sense to try to optimize the distribution of the gaslift if the amount
-    //   of gaslift is constrained either by the maximum allowed gaslift or total gas
-    //   or the group is under individual control
     if ((prod_control != Group::ProductionCMode::NONE) && (prod_control != Group::ProductionCMode::FLD))
     {
-        displayDebugMessage_("optimizing", group_name);
+        if (this->debug) {
+            const std::string msg = fmt::format("optimizing (control = {})", Group::ProductionCMode2String(prod_control));
+            displayDebugMessage_(msg, group_name);
+        }
         auto wells = getGroupGliftWells_(group);
         std::vector<GradPair> inc_grads;
         std::vector<GradPair> dec_grads;
@@ -420,7 +422,10 @@ optimizeGroup_(const Group &group)
         removeSurplusALQ_(group, inc_grads, dec_grads);
     }
     else {
-        displayDebugMessage_("skipping", group_name);
+        if (this->debug) {
+            const std::string msg = fmt::format("skipping (control = {})", Group::ProductionCMode2String(prod_control));
+            displayDebugMessage_(msg, group_name);
+        }
     }
 }
 
@@ -441,12 +446,12 @@ optimizeGroupsRecursive_(const Group &group)
 void
 GasLiftStage2::
 recalculateGradientAndUpdateData_(
-    GradPairItr &grad_itr, const std::string& gr_name_dont_limit, bool increase,
+        GradPairItr &grad_itr, const std::string& gr_name_dont_limit, bool increase,
 
-    //incremental and decremental gradients, if 'grads' are incremental, then
+        //incremental and decremental gradients, if 'grads' are incremental, then
-    // 'other_grads' are decremental, or conversely, if 'grads' are decremental, then
+        // 'other_grads' are decremental, or conversely, if 'grads' are decremental, then
-    // 'other_grads' are incremental
+        // 'other_grads' are incremental
-    std::vector<GradPair> &grads,  std::vector<GradPair> &other_grads)
+        std::vector<GradPair> &grads,  std::vector<GradPair> &other_grads)
 {
     // NOTE: We make a copy of the name string instead of taking a reference
     //   since we may have to erase grad_itr (in the "else" condition below)
@@ -457,26 +462,44 @@ recalculateGradientAndUpdateData_(
     // the grads and other grads are synchronized later
     if(this->stage1_wells_.count(name) > 0) {
         GasLiftSingleWell &gs_well = *(this->stage1_wells_.at(name).get());
-        auto grad = calcIncOrDecGrad_(name, gs_well, gr_name_dont_limit, increase);
+        {
-        if (grad) {
+            auto grad = calcIncOrDecGrad_(name, gs_well, gr_name_dont_limit, increase);
-            grad_itr->second = grad->grad;
+            if (grad) {
-            old_grad = updateGrad_(name, *grad, increase);
+                grad_itr->second = grad->grad;
+                old_grad = updateGrad_(name, *grad, increase);
+            }
+            else {
+                grads.erase(grad_itr); // NOTE: this invalidates grad_itr
+                old_grad = deleteGrad_(name, increase);
+            }
         }
-        else {
+        // If the old incremental/decremental gradient was defined, it becomes the new
-            grads.erase(grad_itr); // NOTE: this invalidates grad_itr
+        //   decremental/incremental gradient
-            old_grad = deleteGrad_(name, increase);
+        if (old_grad) {
+            // NOTE: Either creates a new item or reassigns
+            // The old incremental gradient becomes the new decremental gradient
+            //   or the old decremental gradient becomes the new incremental gradient
+            // NOTE: The gradient itself (old_grad.grad) will be equal to the new decremental
+            //   gradient, but the other fields in old_grad cannot be used as they refer to
+            //   the incremental gradient. E.g. old_grad.alq is the alq after an increase in alq,
+            //   not a decrease, so we need to recalculate the decremental values here..
+            auto grad = calcIncOrDecGrad_(name, gs_well, gr_name_dont_limit, !increase);
+            if (grad) {
+                updateGrad_(name, *grad, !increase);
+                // NOTE: This may invalidate any iterator into 'other_grads' since
+                //   updateGradVector_() will do a push_back() if 'name' is not found..
+                updateGradVector_(name, other_grads, grad->grad);
+            }
+            else {
+                for (auto it = other_grads.begin(); it != other_grads.end(); it++) {
+                    if (it->first == name) {
+                        other_grads.erase(it);
+                        deleteGrad_(name, !increase);
+                    }
+                }
+            }
         }
     }
-
-    if (old_grad) {
-        // NOTE: Either creates a new item or reassigns
-        // The old incremental gradient becomes the new decremental gradient
-        //   or the old decremental gradient becomes the new incremental gradient
-        updateGrad_(name, *old_grad, !increase);
-        // NOTE: This may invalidate any iterator into 'other_grads' since
-        //   updateGradVector_() will do a push_back() if 'name' is not found..
-        updateGradVector_(name, other_grads, old_grad->grad);
-    }
 }
 
 // redistributeALQ_() :
@@ -608,7 +631,7 @@ removeSurplusALQ_(const Group &group,
     if (this->debug) {
         std::string max_glift_str = "unlimited";
         if (max_glift) max_glift_str = fmt::format("{}", *max_glift);
-        const std::string msg = fmt::format("Starting iteration for group: {}. "
+        const std::string msg = fmt::format("Starting remove surplus iteration for group: {}. "
             "oil_rate = {}, oil_target = {}, gas_rate = {}, gas_target = {}, "
             "water_rate = {}, liquid_target = {}, alq = {}, max_alq = {}",
             group.name(), oil_rate, controls.oil_target,
@@ -798,7 +821,7 @@ void
 GasLiftStage2::OptimizeState::
 debugShowIterationInfo()
 {
-    const std::string msg = fmt::format("iteration {}", this->it);
+    const std::string msg = fmt::format("redistribute ALQ iteration {}", this->it);
     displayDebugMessage_(msg);
 }
 

opm/simulators/wells/GasLiftStage2.hpp

@@ -80,7 +80,7 @@ protected:
         GradMap& grad_map, const std::string& well_name, bool add);
     std::optional<GradInfo> calcIncOrDecGrad_(
         const std::string name, const GasLiftSingleWell& gs_well, const std::string& gr_name_dont_limit, bool increase);
-    bool checkRateAlreadyLimited_(GasLiftWellState& state, bool increase);
+    bool checkRateAlreadyLimited_(const std::string& well_name, GasLiftWellState& state, bool increase);
     GradInfo deleteDecGradItem_(const std::string& name);
     GradInfo deleteIncGradItem_(const std::string& name);
     GradInfo deleteGrad_(const std::string& name, bool increase);