From c2916bc7d9011e847474493c50ac3da09a85f03a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Thu, 2 Mar 2023 08:23:15 +0100 Subject: [PATCH 1/2] Fix bug when recalculating gaslift gradient Fixes a bug in recalculateGradientAndUpdateData_() if gaslift stage2. Also adds some more debugging output. --- opm/simulators/wells/GasLiftSingleWell.hpp | 2 +- .../wells/GasLiftSingleWellGeneric.cpp | 34 +++-- .../wells/GasLiftSingleWellGeneric.hpp | 10 +- .../wells/GasLiftSingleWell_impl.hpp | 15 +- opm/simulators/wells/GasLiftStage2.cpp | 141 +++++++++++------- opm/simulators/wells/GasLiftStage2.hpp | 2 +- 6 files changed, 118 insertions(+), 86 deletions(-) diff --git a/opm/simulators/wells/GasLiftSingleWell.hpp b/opm/simulators/wells/GasLiftSingleWell.hpp index d95275a76..11c2d92f2 100644 --- a/opm/simulators/wells/GasLiftSingleWell.hpp +++ b/opm/simulators/wells/GasLiftSingleWell.hpp @@ -55,7 +55,7 @@ namespace Opm const WellInterfaceGeneric &getWell() const override { return well_; } private: - std::optional computeBhpAtThpLimit_(double alq) const override; + std::optional 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); diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp index 1a1f9de89..cb425e42e 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp @@ -85,8 +85,14 @@ GasLiftSingleWellGeneric::GasLiftSingleWellGeneric(DeferredLogger& deferred_logg ****************************************/ // NOTE: Used from GasLiftStage2 std::optional -GasLiftSingleWellGeneric::calcIncOrDecGradient( - double oil_rate, double gas_rate, double water_rate, double alq, const std::string& gr_name_dont_limit, bool increase) const +GasLiftSingleWellGeneric::calcIncOrDecGradient( double oil_rate, + 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); - // 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); + auto rates = computeWellRates_(new_bhp, bhp_is_limited, debug_output); 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 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}; @@ -1161,8 +1163,8 @@ GasLiftSingleWellGeneric::runOptimizeLoop_(bool increase) // if (this->debug) debugShowBhpAlqTable_(); if (this->debug) debugShowAlqIncreaseDecreaseCounts_(); - if (this->debug) - debugShowTargets_(); + //if (this->debug) + // debugShowTargets_(); bool success = false; // did we succeed to increase alq? bool alq_is_limited = false; LimitedRates new_rates = *rates; @@ -1528,10 +1530,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 @@ -1590,7 +1592,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); diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp index c15a0d66d..3d3ea73fb 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp @@ -95,9 +95,11 @@ public: std::optional 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 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 computeBhpAtThpLimit_(double alq) const = 0; + virtual std::optional computeBhpAtThpLimit_(double alq, bool debug_output = true) const = 0; std::pair,double> computeConvergedBhpAtThpLimitByMaybeIncreasingALQ_() const; std::pair,double> computeInitialWellRates_() const; std::optional 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 computeWellRatesWithALQ_(double alq) const; void debugCheckNegativeGradient_(double grad, double alq, double new_alq, double oil_rate, double new_oil_rate, diff --git a/opm/simulators/wells/GasLiftSingleWell_impl.hpp b/opm/simulators/wells/GasLiftSingleWell_impl.hpp index 2267becc3..0f0ec8397 100644 --- a/opm/simulators/wells/GasLiftSingleWell_impl.hpp +++ b/opm/simulators/wells/GasLiftSingleWell_impl.hpp @@ -128,7 +128,7 @@ computeWellRates_( double bhp, bool bhp_is_limited, bool debug_output ) const template std::optional GasLiftSingleWell:: -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( - "Computed bhp ({}) from thp limit is below bhp limit ({}), (ALQ = {})." - " Using bhp limit instead", - *bhp_at_thp_limit, this->controls_.bhp_limit, alq); - displayDebugMessage_(msg); + if (debug_output && this->debug) { + const std::string msg = fmt::format( + "Computed bhp ({}) from thp limit is below bhp limit ({}), (ALQ = {})." + " Using bhp limit instead", + *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); diff --git a/opm/simulators/wells/GasLiftStage2.cpp b/opm/simulators/wells/GasLiftStage2.cpp index d4f15f69f..941dd574a 100644 --- a/opm/simulators/wells/GasLiftStage2.cpp +++ b/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)) { - /* - const std::string msg = fmt::format( - "well {} : not able to obtain {} gradient", - well_name, - (increase ? "incremental" : "decremental") - ); - displayDebugMessage_(msg); - */ + if (checkRateAlreadyLimited_(well_name, state, increase)) { 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); - if (grad) { - const std::string msg = fmt::format( - "well {} : adding {} gradient = {}", - well_name, - (increase ? "incremental" : "decremental"), - grad->grad - ); - displayDebugMessage_(msg); + auto grad = gs_well.calcIncOrDecGradient( + oil_rate, gas_rate, state.waterRate(), alq, gr_name_dont_limit, increase, /*debug_output=*/false); + if (this->debug) { + if (grad) { + const std::string msg = fmt::format( + "well {} : alq = {} : adding {} gradient = {}", + well_name, + alq, + (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.gasIsLimited() ? "gas" : "alq"))); + (state.oilIsLimited() ? "oil" : (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 inc_grads; std::vector 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 - // 'other_grads' are decremental, or conversely, if 'grads' are decremental, then - // 'other_grads' are incremental - std::vector &grads, std::vector &other_grads) + //incremental and decremental gradients, if 'grads' are incremental, then + // 'other_grads' are decremental, or conversely, if 'grads' are decremental, then + // 'other_grads' are incremental + std::vector &grads, std::vector &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,46 @@ 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) { - grad_itr->second = grad->grad; - old_grad = updateGrad_(name, *grad, increase); + { + auto grad = calcIncOrDecGrad_(name, gs_well, gr_name_dont_limit, increase); + if (grad) { + 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 { - grads.erase(grad_itr); // NOTE: this invalidates grad_itr - old_grad = deleteGrad_(name, increase); + // If the old incremental/decremental gradient was defined, it becomes the new + // decremental/incremental gradient + 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); + //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, grad->grad); + //updateGradVector_(name, other_grads, old_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 +633,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 +823,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); } diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp index ae4c022ae..d46c3c15c 100644 --- a/opm/simulators/wells/GasLiftStage2.hpp +++ b/opm/simulators/wells/GasLiftStage2.hpp @@ -80,7 +80,7 @@ protected: GradMap& grad_map, const std::string& well_name, bool add); std::optional 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); From 973a8830435724bfc7915051f750f5d9d3502b87 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Fri, 31 Mar 2023 09:26:28 +0200 Subject: [PATCH 2/2] Some minor changes as requested in review. --- opm/simulators/wells/GasLiftSingleWellGeneric.cpp | 4 ++-- opm/simulators/wells/GasLiftStage2.cpp | 2 -- 2 files changed, 2 insertions(+), 4 deletions(-) diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp index cb425e42e..dbc508151 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp @@ -1163,8 +1163,8 @@ GasLiftSingleWellGeneric::runOptimizeLoop_(bool increase) // if (this->debug) debugShowBhpAlqTable_(); if (this->debug) debugShowAlqIncreaseDecreaseCounts_(); - //if (this->debug) - // debugShowTargets_(); + if (this->debug) + debugShowTargets_(); bool success = false; // did we succeed to increase alq? bool alq_is_limited = false; LimitedRates new_rates = *rates; diff --git a/opm/simulators/wells/GasLiftStage2.cpp b/opm/simulators/wells/GasLiftStage2.cpp index 941dd574a..d1150c7b8 100644 --- a/opm/simulators/wells/GasLiftStage2.cpp +++ b/opm/simulators/wells/GasLiftStage2.cpp @@ -486,11 +486,9 @@ recalculateGradientAndUpdateData_( auto grad = calcIncOrDecGrad_(name, gs_well, gr_name_dont_limit, !increase); if (grad) { updateGrad_(name, *grad, !increase); - //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, grad->grad); - //updateGradVector_(name, other_grads, old_grad->grad); } else { for (auto it = other_grads.begin(); it != other_grads.end(); it++) {