Merge pull request #5694 from totto82/fix_gaslift_PR1

Fix gaslift 1

opm/simulators/wells/GasLiftGroupInfo.cpp

@@ -532,23 +532,30 @@ getProducerWellRates_(const Well* well, int well_index)
 
     const auto controls = well->productionControls(this->summary_state_);
     Scalar oil_rate = oil_pot;
-    if (controls.hasControl(Well::ProducerCMode::ORAT)) {
+    Scalar water_rate = water_pot;
-        oil_rate = std::min(static_cast<Scalar>(controls.oil_rate), oil_rate);
-    }
     Scalar gas_rate = gas_pot;
+
+    if (controls.hasControl(Well::ProducerCMode::ORAT) && oil_rate > static_cast<Scalar>(controls.oil_rate)) {
+        water_rate *= (static_cast<Scalar>(controls.oil_rate) / oil_rate); 
+        oil_rate = static_cast<Scalar>(controls.oil_rate);
+    }
+    
     if (controls.hasControl(Well::ProducerCMode::GRAT)) {
         gas_rate = std::min(static_cast<Scalar>(controls.gas_rate), gas_rate);
     }
-    Scalar water_rate = water_pot;
+    
-    if (controls.hasControl(Well::ProducerCMode::WRAT)) {
+    if (controls.hasControl(Well::ProducerCMode::WRAT) && water_rate > static_cast<Scalar>(controls.water_rate)) {
-        water_rate = std::min(static_cast<Scalar>(controls.water_rate), water_rate);
+        oil_rate *= (static_cast<Scalar>(controls.water_rate) / water_rate); 
+        water_rate = static_cast<Scalar>(controls.water_rate);
     }
     if (controls.hasControl(Well::ProducerCMode::LRAT)) {
         Scalar liquid_rate = oil_rate + water_rate;
-        Scalar liquid_rate_lim = std::min(static_cast<Scalar>(controls.liquid_rate), liquid_rate);
+        Scalar liquid_rate_lim = static_cast<Scalar>(controls.liquid_rate);
+        if (liquid_rate > liquid_rate_lim) {
             water_rate = water_rate / liquid_rate * liquid_rate_lim;
             oil_rate = oil_rate / liquid_rate * liquid_rate_lim;
         }
+    }
 
     return {oil_rate, gas_rate, water_rate, oil_pot, gas_pot, water_pot};
 }
@@ -639,19 +646,6 @@ initializeGroupRatesRecursive_(const Group& group)
     }
     if (oil_target || liquid_target || water_target || gas_target || max_total_gas || max_alq) {
         updateGroupIdxMap_(group.name());
-        if (oil_target)
-            oil_rate = std::min(oil_rate, *oil_target);
-        if (gas_target)
-            gas_rate = std::min(gas_rate, *gas_target);
-        if (water_target)
-            water_rate = std::min(water_rate, *water_target);
-        if (liquid_target) {
-            Scalar liquid_rate = oil_rate + water_rate;
-            Scalar liquid_rate_limited = std::min(liquid_rate, *liquid_target);
-            oil_rate = oil_rate / liquid_rate * liquid_rate_limited;
-            water_rate = water_rate / liquid_rate * liquid_rate_limited;
-        }
-
         this->group_rate_map_.try_emplace(group.name(),
             oil_rate, gas_rate, water_rate, alq,
             oil_potential, gas_potential, water_potential,
@@ -660,6 +654,25 @@ initializeGroupRatesRecursive_(const Group& group)
             debugDisplayUpdatedGroupRates(
                 group.name(), oil_rate, gas_rate, water_rate, alq);
         }
+
+        if (oil_target && oil_rate > *oil_target)  {
+            water_rate *= (*oil_target/oil_rate);
+            oil_rate = *oil_target;
+        }
+        if (gas_target && gas_rate > *gas_target)
+            gas_rate = *gas_target;
+        if (water_target && water_rate > *water_target) {
+            oil_rate *= (*water_target/water_rate); 
+            water_rate = *water_target;
+        }
+        if (liquid_target) {
+            Scalar liquid_rate = oil_rate + water_rate;
+            Scalar liquid_rate_limited = *liquid_target;
+            if (liquid_rate > liquid_rate_limited) {
+                oil_rate = oil_rate / liquid_rate * liquid_rate_limited;
+                water_rate = water_rate / liquid_rate * liquid_rate_limited;
+            }
+        }
     }
     return std::make_tuple(oil_rate, gas_rate, water_rate, oil_potential, gas_potential, water_potential, alq);
 }
@@ -687,7 +700,9 @@ initializeWell2GroupMapRecursive_(const Group& group,
             // TODO: can the same well be memember of two different groups
             //  (on the same recursion level) ?
             assert(this->well_group_map_.count(well_name) == 0);
-            if (checkDoGasLiftOptimization_(well_name)) {
+            bool checkDoGasLift = checkDoGasLiftOptimization_(well_name);
+            checkDoGasLift = this->comm_.max(checkDoGasLift);
+            if (checkDoGasLift) {
                 const auto &well = this->schedule_.getWell(
                     well_name, this->report_step_idx_);
                 Scalar wfac = well.getEfficiencyFactor();

opm/simulators/wells/GasLiftSingleWellGeneric.cpp

@@ -204,7 +204,7 @@ addOrSubtractAlqIncrement_(Scalar alq, bool increase) const
                 limited = true;
             }
         } else {
-            if (alq < 0) {
+            if (alq < 1e-12) {
                 alq = 0.0;
                 limited = true;
             }
@@ -261,7 +261,7 @@ checkGroupTargetsViolated(const BasicRates& rates,
             auto target_opt = this->group_info_.getTarget(rate_type, group_name);
             if (target_opt) {
                 auto delta_rate = new_rates[rate_type] - rates[rate_type];
-                auto new_group_rate = this->group_info_.getPotential(rate_type, group_name) + efficiency * delta_rate;
+                auto new_group_rate = this->group_info_.getRate(rate_type, group_name) + efficiency * delta_rate;
                 if (new_group_rate > *target_opt) {
                     if (this->debug) {
                         const std::string msg
@@ -666,6 +666,14 @@ getRateWithLimit_(Rate rate_type, const BasicRates& rates) const
     //   for why the rate was limited.
     std::optional<Rate> target_type;
 
+    // we also need to limit the other rate (currently only for water and oil and not gas)
+    Scalar rate2 = 0.0;
+    if (rate_type == Rate::oil) {
+        rate2 = getRate_(Rate::water, rates);
+    } else if (rate_type == Rate::water) {
+        rate2 = getRate_(Rate::oil, rates);
+    }
+
     if (hasProductionControl_(rate_type)) {
         auto target = getProductionTarget_(rate_type);
         if (new_rate > target) {
@@ -675,17 +683,47 @@ getRateWithLimit_(Rate rate_type, const BasicRates& rates) const
                                                 new_rate,
                                                 target);
             displayDebugMessage_(msg);
+            rate2 *= target/new_rate;
             new_rate = target;
             target_type = rate_type; 
         }
     }
-    if (((rate_type == Rate::oil) || (rate_type == Rate::water)) && hasProductionControl_(Rate::liquid)) {
+    if (((rate_type == Rate::oil) || (rate_type == Rate::water))) {
-        Scalar rate2;
         if (rate_type == Rate::oil) {
-            rate2 = getRate_(Rate::water, rates);
+            if(hasProductionControl_(Rate::water)) {
-        } else {
+                auto water_target = getProductionTarget_(Rate::water);
-            rate2 = getRate_(Rate::oil, rates);
+                if (rate2 > water_target) {
+                    new_rate *= (water_target / rate2);
+                    target_type = Rate::water;
+                    rate2 = water_target;
+                    const std::string msg = fmt::format("limiting {} rate to {} due to WRAT target: "
+                                                        "computed WRAT: {}, target WRAT: {}",
+                                                        GasLiftGroupInfo<Scalar>::rateToString(rate_type),
+                                                        new_rate,
+                                                        rate2,
+                                                        water_target);
+                    displayDebugMessage_(msg);
                 }
+            }
+        } else {
+            if(hasProductionControl_(Rate::oil)) {
+                auto oil_target = getProductionTarget_(Rate::oil);
+                if (rate2 > oil_target) {
+                    new_rate *= (oil_target / rate2);
+                    target_type = Rate::oil;
+                    rate2 = oil_target;
+                    const std::string msg = fmt::format("limiting {} rate to {} due to ORAT target: "
+                                                        "computed ORAT: {}, target ORAT: {}",
+                                                        GasLiftGroupInfo<Scalar>::rateToString(rate_type),
+                                                        new_rate,
+                                                        rate2,
+                                                        oil_target);
+                    displayDebugMessage_(msg);
+                }
+            }
+        }
+
+        if(hasProductionControl_(Rate::liquid)) {
             // Note: Since "new_rate" was first updated for ORAT or WRAT, see first "if"
             //   statement in the method, the rate is limited due to LRAT only if
             //   it becomes less than the rate limited by a WRAT or ORAT target..
@@ -713,6 +751,7 @@ getRateWithLimit_(Rate rate_type, const BasicRates& rates) const
                 displayDebugMessage_(msg);
             }
         }
+    }
     // TODO: Also check RESV target?
     return {new_rate, target_type};
 }
@@ -768,25 +807,7 @@ getLiquidRateWithGroupLimit_(const Scalar new_oil_rate,
         = getRateWithGroupLimit_(Rate::liquid, new_liquid_rate, liquid_rate, gr_name_dont_limit);
     bool limited = group_name != nullptr;
     if (limited) {
-        // the oil, gas, and water cases can be handled directly by
+        Scalar oil_fraction = oil_rate / liquid_rate;
-        //  getRateWithGroupLimit_() above. However, for the liquid case
-        //  we must do some postprocessing. I chose to include it here
-        //  instead of cluttering up getRateWithGroupLimit_() with this
-        //  special case.
-        Scalar delta_water = new_water_rate - water_rate;
-        Scalar delta_oil = new_oil_rate - oil_rate;
-
-        Scalar gr_water_rate = this->group_info_.waterRate(*group_name);
-        Scalar gr_oil_rate = this->group_info_.oilRate(*group_name);
-
-        // NOTE: these rates are too large according to the limited liquid rate
-        //  but it does not matter since we are only using them to calculate
-        //  the fraction of the liquid corresponding to the oil phase
-        Scalar new_gr_water_rate = gr_water_rate + efficiency * delta_water;
-        Scalar new_gr_oil_rate = gr_oil_rate + efficiency * delta_oil;
-        Scalar new_gr_liquid_rate = new_gr_water_rate + new_gr_oil_rate;
-
-        Scalar oil_fraction = new_gr_oil_rate / new_gr_liquid_rate;
         Scalar delta_liquid = liquid_rate_limited - liquid_rate;
         auto limited_oil_rate = oil_rate + oil_fraction * delta_liquid;
         auto limited_water_rate = water_rate + (1.0 - oil_fraction) * delta_liquid;
@@ -880,6 +901,7 @@ getInitialRatesWithLimit_() const
         auto temp_rates = getLimitedRatesFromRates_(*rates);
         BasicRates old_rates = getWellStateRates_();
         limited_rates = updateRatesToGroupLimits_(old_rates, temp_rates);
+
         initial_alq = alq;
     }
     return {limited_rates, initial_alq};
@@ -891,9 +913,10 @@ GasLiftSingleWellGeneric<Scalar>::
 getLimitedRatesFromRates_(const BasicRates& rates) const
 {
     auto [oil_rate, oil_limiting_target] = getOilRateWithLimit2_(rates);
+    bool oil_is_limited = oil_limiting_target.has_value();
     auto [gas_rate, gas_is_limited] = getGasRateWithLimit_(rates);
     auto [water_rate, water_limiting_target] = getWaterRateWithLimit2_(rates);
-    bool oil_is_limited = oil_limiting_target.has_value();
+   
     bool water_is_limited = water_limiting_target.has_value();
     return LimitedRates {oil_rate,
                          gas_rate,
@@ -1071,8 +1094,6 @@ maybeAdjustALQbeforeOptimizeLoop_(const LimitedRates& orig_rates,
             const std::string msg = fmt::format("adjusted ALQ to: {}", alq);
             displayDebugMessage_(msg);
         }
-        Scalar delta_alq = alq - orig_alq;
-        updateGroupRates_(orig_rates, rates, delta_alq);
     }
     return {rates, alq};
 }
@@ -1465,13 +1486,16 @@ updateRatesToGroupLimits_(const BasicRates& old_rates,
 {
     LimitedRates new_rates = rates;
     auto [new_oil_rate, oil_is_limited] = getOilRateWithGroupLimit_(new_rates.oil, old_rates.oil, gr_name);
+    auto mod_water_rate = new_rates.water;
     if (oil_is_limited) {
         new_rates.oil_limiting_target = Rate::oil;
+        mod_water_rate *=  (new_oil_rate / new_rates.oil);
     }
     auto [new_gas_rate, gas_is_limited] = getGasRateWithGroupLimit_(new_rates.gas, old_rates.gas, gr_name);
-    auto [new_water_rate, water_is_limited] = getWaterRateWithGroupLimit_(new_rates.water, old_rates.water, gr_name);
+    auto [new_water_rate, water_is_limited] = getWaterRateWithGroupLimit_(mod_water_rate, old_rates.water, gr_name);
     if (water_is_limited) {
         new_rates.water_limiting_target = Rate::water;
+        new_oil_rate *=  (new_water_rate / new_rates.water);
     }
     auto [new_oil_rate2, new_water_rate2, oil_is_limited2, water_is_limited2]
         = getLiquidRateWithGroupLimit_(new_oil_rate, old_rates.oil, new_water_rate, old_rates.water, gr_name);

opm/simulators/wells/GasLiftStage2.cpp

@@ -117,11 +117,6 @@ addOrRemoveALQincrement_(GradMap &grad_map,
             well_name, (add ? "adding" : "subtracting"), old_alq, new_alq);
         this->displayDebugMessage_(msg);
     }
-    state.update(gi.new_oil_rate, gi.oil_is_limited,
-                 gi.new_gas_rate, gi.gas_is_limited,
-                 gi.alq, gi.alq_is_limited,
-                 gi.new_water_rate, gi.water_is_limited, add);
-
     this->well_state_.setALQ(well_name, gi.alq);
     const auto& pu = this->well_state_.phaseUsage();
     std::vector<Scalar> well_pot(pu.num_phases, 0.0);
@@ -426,7 +421,7 @@ optimizeGroup_(const Group& group)
 {
     const auto& group_name = group.name();
     const auto prod_control = this->group_state_.production_control(group_name);
-    if (this->glo_.has_group(group.name()) || ((prod_control != Group::ProductionCMode::NONE) && (prod_control != Group::ProductionCMode::FLD)))
+    if ((prod_control != Group::ProductionCMode::NONE) && (prod_control != Group::ProductionCMode::FLD))
     {
         if (this->debug) {
             const std::string msg = fmt::format("optimizing (control = {})", Group::ProductionCMode2String(prod_control));
@@ -436,7 +431,7 @@ optimizeGroup_(const Group& group)
         std::vector<GradPair> inc_grads;
         std::vector<GradPair> dec_grads;
         redistributeALQ_(wells, group, inc_grads, dec_grads);
-        removeSurplusALQ_(group, inc_grads, dec_grads);
+        removeSurplusALQ_(group, dec_grads);
     }
     else {
         if (this->debug) {
@@ -493,7 +488,7 @@ recalculateGradientAndUpdateData_(GradPairItr& grad_itr,
         }
         // If the old incremental/decremental gradient was defined, it becomes the new
         //   decremental/incremental gradient
-        if (old_grad) {
+        if (old_grad && !other_grads.empty()) {
             // 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
@@ -634,13 +629,13 @@ redistributeALQ_(std::vector<GasLiftSingleWell*>& wells,
 template<class Scalar>
 void GasLiftStage2<Scalar>::
 removeSurplusALQ_(const Group& group,
-                  std::vector<GradPair>& inc_grads,
                   std::vector<GradPair>& dec_grads)
 {
     if (dec_grads.empty()) {
         displayDebugMessage_("no wells to remove ALQ from. Skipping");
         return;
     }
+    std::vector<GradPair> empty_vector;
     assert(!dec_grads.empty());
     const auto max_glift = getGroupMaxALQ_(group);
     const auto max_totalgas = getGroupMaxTotalGas_(group);
@@ -697,14 +692,14 @@ removeSurplusALQ_(const Group& group,
             state.addOrRemoveALQincrement( this->dec_grads_, well_name, /*add=*/false);
             // We pass the current group rate in order to avoid limiting the rates
             // and gaslift based on the current group limits. In other words we want to reduce
-            // the gasslift as much as possible as long as we are able to produce the group
+            // the gaslift as much as possible as long as we are able to produce the group
             // targets
+            // Note: empty_vector passed to avoid adding gradients to inc_grads
             recalculateGradientAndUpdateData_(
-                        dec_grad_itr, group.name(), /*increase=*/false, dec_grads, inc_grads);
+                        dec_grad_itr, group.name(), /*increase=*/false, dec_grads, empty_vector);
 
-            // The dec_grads and inc_grads needs to be syncronized across ranks
+            // The dec_grads needs to be syncronized across ranks
             mpiSyncGlobalGradVector_(dec_grads);
-            mpiSyncGlobalGradVector_(inc_grads);
             // NOTE: recalculateGradientAndUpdateData_() will remove the current gradient
             //   from dec_grads if it cannot calculate a new decremental gradient.
             //   This will invalidate dec_grad_itr and well_name
@@ -798,6 +793,59 @@ updateGradVector_(const std::string& name,
     //   later in getEcoGradients()
 }
 
+
+template<class Scalar>
+void
+GasLiftStage2<Scalar>::
+updateGroupInfo(const std::string& well_name, bool add)
+{
+    const auto delta = computeDelta(well_name, add);
+    const auto& [delta_oil, delta_gas, delta_water, delta_alq] = delta;
+    if (this->group_info_.hasWell(well_name)) {
+        const auto& pairs = this->group_info_.getWellGroups(well_name);
+       for (const auto& [group_name, efficiency] : pairs) {
+            this->group_info_.update(group_name,
+                                     efficiency * delta_oil,
+                                     efficiency * delta_gas,
+                                     efficiency * delta_water,
+                                     efficiency * delta_alq);
+        }
+    }
+}
+
+template<class Scalar>
+std::array<Scalar, 4>
+GasLiftStage2<Scalar>::
+computeDelta(const std::string& well_name, bool add)
+{
+    std::array<Scalar, 4> delta = {0.0, 0.0, 0.0, 0.0};
+    // compute the delta on wells on own rank
+    if (this->well_state_map_.count(well_name) > 0) {
+        const GradInfo& gi = add? this->inc_grads_.at(well_name) : this->dec_grads_.at(well_name);
+        GasLiftWellState<Scalar>& state = *(this->well_state_map_.at(well_name).get());
+        GasLiftSingleWell& gs_well = *(this->stage1_wells_.at(well_name).get());
+        const WellInterfaceGeneric<Scalar>& well = gs_well.getWell();
+        // only get deltas for wells owned by this rank
+        if (this->well_state_.wellIsOwned(well.indexOfWell(), well_name)) {
+            const auto& well_ecl = well.wellEcl();
+            Scalar factor = well_ecl.getEfficiencyFactor();
+            auto& [delta_oil, delta_gas, delta_water, delta_alq] = delta;
+            delta_oil = factor * (gi.new_oil_rate - state.oilRate());
+            delta_gas = factor * (gi.new_gas_rate - state.gasRate());
+            delta_water = factor * (gi.new_water_rate - state.waterRate());
+            delta_alq = factor * (gi.alq - state.alq());
+        }
+        state.update(gi.new_oil_rate, gi.oil_is_limited,
+                gi.new_gas_rate, gi.gas_is_limited,
+                gi.alq, gi.alq_is_limited,
+                gi.new_water_rate, gi.water_is_limited, add);
+    }
+
+    // and communicate the results
+    this->comm_.sum(delta.data(), delta.size());
+
+    return delta;
+}
 /***********************************************
  * Public methods declared in OptimizeState
  ***********************************************/
@@ -930,6 +978,9 @@ redistributeALQ(GradPairItr& min_dec_grad,
         this->parent.dec_grads_, /*well_name=*/min_dec_grad->first, /*add=*/false);
     this->parent.addOrRemoveALQincrement_(
         this->parent.inc_grads_, /*well_name=*/max_inc_grad->first, /*add=*/true);
+
+    this->parent.updateGroupInfo(min_dec_grad->first, /*add=*/false);
+    this->parent.updateGroupInfo(max_inc_grad->first, /*add=*/true);
 }
 
 /**********************************************
@@ -966,6 +1017,7 @@ addOrRemoveALQincrement(GradMap& grad_map,
         this->parent.displayDebugMessage_(msg);
     }
     this->parent.addOrRemoveALQincrement_(grad_map, well_name, add);
+    this->parent.updateGroupInfo(well_name, add);
 }
 
 template<class Scalar>

opm/simulators/wells/GasLiftStage2.hpp

@@ -127,7 +127,6 @@ protected:
                           std::vector<GradPair>& dec_grads);
 
     void removeSurplusALQ_(const Group& group,
-                           std::vector<GradPair>& inc_grads,
                            std::vector<GradPair>& dec_grads);
 
     void saveGrad_(GradMap& map, const std::string& name, GradInfo& grad);
@@ -146,6 +145,10 @@ protected:
     void mpiSyncLocalToGlobalGradVector_(const std::vector<GradPair>& grads_local,
                                          std::vector<GradPair>& grads_global) const;
 
+    std::array<Scalar, 4> computeDelta(const std::string& name, bool add);
+    void updateGroupInfo(const std::string& name, bool add);
+
+
     GLiftProdWells& prod_wells_;
     GLiftOptWells& stage1_wells_;
     GasLiftGroupInfo<Scalar>& group_info_;