diff --git a/opm/simulators/wells/GasLiftGroupInfo.cpp b/opm/simulators/wells/GasLiftGroupInfo.cpp index 8222fa3f9..248531186 100644 --- a/opm/simulators/wells/GasLiftGroupInfo.cpp +++ b/opm/simulators/wells/GasLiftGroupInfo.cpp @@ -532,22 +532,29 @@ 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)) { - oil_rate = std::min(static_cast(controls.oil_rate), oil_rate); - } + Scalar water_rate = water_pot; Scalar gas_rate = gas_pot; + + if (controls.hasControl(Well::ProducerCMode::ORAT) && oil_rate > static_cast(controls.oil_rate)) { + water_rate *= (static_cast(controls.oil_rate) / oil_rate); + oil_rate = static_cast(controls.oil_rate); + } + if (controls.hasControl(Well::ProducerCMode::GRAT)) { gas_rate = std::min(static_cast(controls.gas_rate), gas_rate); } - Scalar water_rate = water_pot; - if (controls.hasControl(Well::ProducerCMode::WRAT)) { - water_rate = std::min(static_cast(controls.water_rate), water_rate); + + if (controls.hasControl(Well::ProducerCMode::WRAT) && water_rate > static_cast(controls.water_rate)) { + oil_rate *= (static_cast(controls.water_rate) / water_rate); + water_rate = static_cast(controls.water_rate); } if (controls.hasControl(Well::ProducerCMode::LRAT)) { Scalar liquid_rate = oil_rate + water_rate; - Scalar liquid_rate_lim = std::min(static_cast(controls.liquid_rate), liquid_rate); - water_rate = water_rate / liquid_rate * liquid_rate_lim; - oil_rate = oil_rate / liquid_rate * liquid_rate_lim; + Scalar liquid_rate_lim = static_cast(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(); diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp index e87e5a63c..c77997180 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp +++ b/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 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,42 +683,73 @@ getRateWithLimit_(Rate rate_type, const BasicRates& rates) const new_rate, target); displayDebugMessage_(msg); + rate2 *= target/new_rate; new_rate = target; - target_type = rate_type; + target_type = rate_type; } } - if (((rate_type == Rate::oil) || (rate_type == Rate::water)) && hasProductionControl_(Rate::liquid)) { - Scalar rate2; + if (((rate_type == Rate::oil) || (rate_type == Rate::water))) { if (rate_type == Rate::oil) { - rate2 = getRate_(Rate::water, rates); + if(hasProductionControl_(Rate::water)) { + auto water_target = getProductionTarget_(Rate::water); + 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::rateToString(rate_type), + new_rate, + rate2, + water_target); + displayDebugMessage_(msg); + } + } } else { - rate2 = getRate_(Rate::oil, rates); + 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::rateToString(rate_type), + new_rate, + rate2, + oil_target); + displayDebugMessage_(msg); + } + } } - // 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.. - Scalar liq_rate = new_rate + rate2; - auto liq_target = getProductionTarget_(Rate::liquid); - if (liq_rate > liq_target) { - Scalar fraction = new_rate / liq_rate; - // NOTE: since - // fraction * liq_rate = new_rate, - // we must have - // fraction * liq_target < new_rate - // since - // liq_target < liq_rate - // therefore new_rate will become less than it original was and - // limited = true. - new_rate = fraction * liq_target; - target_type = Rate::liquid; - const std::string msg = fmt::format("limiting {} rate to {} due to LRAT target: " - "computed LRAT: {}, target LRAT: {}", - GasLiftGroupInfo::rateToString(rate_type), - new_rate, - liq_rate, - liq_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.. + Scalar liq_rate = new_rate + rate2; + + auto liq_target = getProductionTarget_(Rate::liquid); + if (liq_rate > liq_target) { + Scalar fraction = new_rate / liq_rate; + // NOTE: since + // fraction * liq_rate = new_rate, + // we must have + // fraction * liq_target < new_rate + // since + // liq_target < liq_rate + // therefore new_rate will become less than it original was and + // limited = true. + new_rate = fraction * liq_target; + target_type = Rate::liquid; + const std::string msg = fmt::format("limiting {} rate to {} due to LRAT target: " + "computed LRAT: {}, target LRAT: {}", + GasLiftGroupInfo::rateToString(rate_type), + new_rate, + liq_rate, + liq_target); + displayDebugMessage_(msg); + } } } // TODO: Also check RESV target? @@ -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 - // 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 oil_fraction = oil_rate / 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:: 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); diff --git a/opm/simulators/wells/GasLiftStage2.cpp b/opm/simulators/wells/GasLiftStage2.cpp index b4f916556..176b2efd9 100644 --- a/opm/simulators/wells/GasLiftStage2.cpp +++ b/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 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 inc_grads; std::vector 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& wells, template void GasLiftStage2:: removeSurplusALQ_(const Group& group, - std::vector& inc_grads, std::vector& dec_grads) { if (dec_grads.empty()) { displayDebugMessage_("no wells to remove ALQ from. Skipping"); return; } + std::vector 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 +void +GasLiftStage2:: +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 +std::array +GasLiftStage2:: +computeDelta(const std::string& well_name, bool add) +{ + std::array 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& state = *(this->well_state_map_.at(well_name).get()); + GasLiftSingleWell& gs_well = *(this->stage1_wells_.at(well_name).get()); + const WellInterfaceGeneric& 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 diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp index e0f5ae3d8..8f28040b9 100644 --- a/opm/simulators/wells/GasLiftStage2.hpp +++ b/opm/simulators/wells/GasLiftStage2.hpp @@ -127,7 +127,6 @@ protected: std::vector& dec_grads); void removeSurplusALQ_(const Group& group, - std::vector& inc_grads, std::vector& dec_grads); void saveGrad_(GradMap& map, const std::string& name, GradInfo& grad); @@ -146,6 +145,10 @@ protected: void mpiSyncLocalToGlobalGradVector_(const std::vector& grads_local, std::vector& grads_global) const; + std::array computeDelta(const std::string& name, bool add); + void updateGroupInfo(const std::string& name, bool add); + + GLiftProdWells& prod_wells_; GLiftOptWells& stage1_wells_; GasLiftGroupInfo& group_info_;