diff --git a/opm/simulators/flow/BlackoilModelEbos.hpp b/opm/simulators/flow/BlackoilModelEbos.hpp index 856416441..85309b774 100644 --- a/opm/simulators/flow/BlackoilModelEbos.hpp +++ b/opm/simulators/flow/BlackoilModelEbos.hpp @@ -416,11 +416,13 @@ namespace Opm { /// \param[in] reservoir_state reservoir state variables /// \param[in, out] well_state well state variables /// \param[in] initial_assembly pass true if this is the first call to assemble() in this timestep - SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface& /* timer */, + SimulatorReportSingle assembleReservoir(const SimulatorTimerInterface &/*timer */, const int iterationIdx) { // -------- Mass balance equations -------- + //auto stepNum = timer.currentStepNum(); ebosSimulator_.model().newtonMethod().setIterationIndex(iterationIdx); + //ebosSimulator_.model().newtonMethod().setCurrentTimeStep(stepNum); ebosSimulator_.problem().beginIteration(); ebosSimulator_.model().linearizer().linearizeDomain(); ebosSimulator_.problem().endIteration(); diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 8ec4adc42..1b36752ce 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -51,10 +51,13 @@ #include #include +#include +#include +#include +#include #include #include #include -#include #include #include #include @@ -103,6 +106,15 @@ namespace Opm { using SparseMatrixAdapter = GetPropType; typedef typename Opm::BaseAuxiliaryModule::NeighborSet NeighborSet; + using GasLiftSingleWell = Opm::GasLiftSingleWell; + using GasLiftStage2 = Opm::GasLiftStage2; + using GLiftWellState = Opm::GasLiftWellState; + using GLiftWellStateMap = + std::map>; + using GLiftOptWells = + std::map>; + using GLiftProdWells = + std::map *>; static const int numEq = Indices::numEq; static const int solventSaturationIdx = Indices::solventSaturationIdx; @@ -346,6 +358,8 @@ namespace Opm { // Check if well equations is converged. ConvergenceReport getWellConvergence(const std::vector& B_avg, const bool checkGroupConvergence = false) const; + const PhaseUsage& phaseUsage() const { return phase_usage_; } + const SimulatorReportSingle& lastReport() const; void addWellContributions(SparseMatrixAdapter& jacobian) const @@ -532,6 +546,14 @@ namespace Opm { void assembleWellEq(const double dt, Opm::DeferredLogger& deferred_logger); + void maybeDoGasLiftOptimize(Opm::DeferredLogger& deferred_logger); + + void gliftDebugShowALQ(Opm::DeferredLogger& deferred_logger); + + void gasLiftOptimizationStage2(Opm::DeferredLogger& deferred_logger, + GLiftProdWells &prod_wells, GLiftOptWells &glift_wells, + GLiftWellStateMap &map); + // some preparation work, mostly related to group control and RESV, // at the beginning of each time step (Not report step) void prepareTimeStep(Opm::DeferredLogger& deferred_logger); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index 416d5b65e..7a0868acb 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -306,8 +306,8 @@ namespace Opm { Opm::DeferredLogger local_deferredLogger; this->resetWellState(); - this->wellState().disableGliftOptimization(); updateAndCommunicateGroupData(); + well_state_.gliftTimeStepInit(); const int reportStepIdx = ebosSimulator_.episodeIndex(); const double simulationTime = ebosSimulator_.time(); std::string exc_msg; @@ -1025,10 +1025,8 @@ namespace Opm { updateWellControls(local_deferredLogger, /* check group controls */ false); } - gliftDebug("assemble() : running assembleWellEq()..", local_deferredLogger); - this->wellState().enableGliftOptimization(); + maybeDoGasLiftOptimize(local_deferredLogger); assembleWellEq(dt, local_deferredLogger); - this->wellState().disableGliftOptimization(); } catch (const std::runtime_error& e) { exc_type = ExceptionType::RUNTIME_ERROR; exc_msg = e.what(); @@ -1047,13 +1045,64 @@ namespace Opm { last_report_.assemble_time_well += perfTimer.stop(); } + maybeDoGasLiftOptimize(Opm::DeferredLogger& deferred_logger) + { + well_state_.enableGliftOptimization(); + GLiftOptWells glift_wells; + GLiftProdWells prod_wells; + GLiftWellStateMap state_map; + // Stage1: Optimize single wells not checking any group limits + for (auto& well : well_container_) { + well->gasLiftOptimizationStage1( + well_state_, ebosSimulator_, deferred_logger, + prod_wells, glift_wells, state_map); + } + gasLiftOptimizationStage2(deferred_logger, prod_wells, glift_wells, state_map); + if (this->glift_debug) gliftDebugShowALQ(deferred_logger); + well_state_.disableGliftOptimization(); + } + + // If a group has any production rate constraints, and/or a limit + // on its total rate of lift gas supply, allocate lift gas + // preferentially to the wells that gain the most benefit from + // it. Lift gas increments are allocated in turn to the well that + // currently has the largest weighted incremental gradient. The + // procedure takes account of any limits on the group production + // rate or lift gas supply applied to any level of group. + template + void + BlackoilWellModel:: + gasLiftOptimizationStage2(Opm::DeferredLogger& deferred_logger, + GLiftProdWells &prod_wells, GLiftOptWells &glift_wells, + GLiftWellStateMap &glift_well_state_map) + { + + GasLiftStage2 glift {*this, ebosSimulator_, deferred_logger, well_state_, + prod_wells, glift_wells, glift_well_state_map}; + glift.runOptimize(); + } + + template + void + BlackoilWellModel:: + gliftDebugShowALQ(Opm::DeferredLogger& deferred_logger) + { + for (auto& well : this->well_container_) { + if (well->isProducer()) { + auto alq = this->well_state_.getALQ(well->name()); + const std::string msg = fmt::format("ALQ_REPORT : {} : {}", + well->name(), alq); + gliftDebug(msg, deferred_logger); + } + } + } + template void BlackoilWellModel:: assembleWellEq(const double dt, Opm::DeferredLogger& deferred_logger) { for (auto& well : well_container_) { - well->maybeDoGasLiftOptimization(this->wellState(), ebosSimulator_, deferred_logger); well->assembleWellEq(ebosSimulator_, dt, this->wellState(), deferred_logger); } } diff --git a/opm/simulators/wells/GasLiftRuntime.hpp b/opm/simulators/wells/GasLiftRuntime.hpp deleted file mode 100644 index 08301bd83..000000000 --- a/opm/simulators/wells/GasLiftRuntime.hpp +++ /dev/null @@ -1,153 +0,0 @@ -/* - Copyright 2020 Equinor ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#ifndef OPM_GASLIFT_RUNTIME_HEADER_INCLUDED -#define OPM_GASLIFT_RUNTIME_HEADER_INCLUDED - -#include -#include -#include -#include -#include -#include -#include -#include -// NOTE: StandardWell.hpp includes ourself (GasLiftRuntime.hpp), so we need -// to forward declare StandardWell for it to be defined in this file. -namespace Opm { - template class StandardWell; -} -#include - -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -namespace Opm -{ - template - class GasLiftRuntime { - using Simulator = GetPropType; - using WellState = WellStateFullyImplicitBlackoil; - using StdWell = Opm::StandardWell; - // TODO: same definition with WellInterface, and - // WellStateFullyImplicitBlackoil eventually they should go - // to a common header file. - static const int Water = BlackoilPhases::Aqua; - static const int Oil = BlackoilPhases::Liquid; - static const int Gas = BlackoilPhases::Vapour; - struct OptimizeState; - public: - GasLiftRuntime( - const StdWell &std_well, - const Simulator &ebos_simulator, - const SummaryState &summary_state, - DeferredLogger &deferred_logger, - WellState &well_state, - const Well::ProductionControls &controls - ); - void runOptimize(); - private: - void computeInitialWellRates_(); - void computeWellRates_(double bhp, std::vector &potentials); - void debugShowIterationInfo_(OptimizeState &state, double alq); - void debugShowStartIteration_(double alq, bool increase); - void displayDebugMessage_(const std::string &msg); - void displayWarning_(); - void displayWarning_(std::string warning); - bool getGasRateWithLimit_(double& new_rate, const std::vector &potentials); - bool getOilRateWithLimit_(double& new_rate, const std::vector &potentials); - void logSuccess_(); - bool runOptimizeLoop_(bool increase); - void setAlqMaxRate_(const GasLiftOpt::Well &well); - void setAlqMinRate_(const GasLiftOpt::Well &well); - bool tryIncreaseLiftGas_(); - bool tryDecreaseLiftGas_(); - void updateWellStateAlqFixedValue_(const GasLiftOpt::Well &well); - bool useFixedAlq_(const GasLiftOpt::Well &well); - void warnMaxIterationsExceeded_(); - - const Well::ProductionControls &controls_; - DeferredLogger &deferred_logger_; - const Simulator &ebos_simulator_; - std::vector potentials_; - const StdWell &std_well_; - const SummaryState &summary_state_; - WellState &well_state_; - std::string well_name_; - bool debug; // extra debug output - - double alpha_w_; - double alpha_g_; - double eco_grad_; - int gas_pos_; - bool has_run_init_ = false; - double increment_; - double max_alq_; - int max_iterations_; - double min_alq_; - double new_alq_; - int oil_pos_; - bool optimize_; - double orig_alq_; - int water_pos_; - - struct OptimizeState { - OptimizeState( GasLiftRuntime &parent_, bool increase_ ) : - parent(parent_), - increase(increase_), - it(0), - stop_iteration(false), - bhp(-1) - {} - - GasLiftRuntime &parent; - bool increase; - int it; - bool stop_iteration; - double bhp; - - double addOrSubtractAlqIncrement(double alq); - double calcGradient(double oil_rate, double new_oil_rate, - double gas_rate, double new_gas_rate); - bool checkAlqOutsideLimits(double alq, double oil_rate); - bool checkEcoGradient(double gradient); - bool checkOilRateExceedsTarget(double oil_rate); - bool checkRate(double rate, double limit, const std::string rate_str); - bool checkWellRatesViolated(std::vector &potentials); - bool computeBhpAtThpLimit(double alq); - double getBhpWithLimit(); - void warn_(std::string msg) {parent.displayWarning_(msg);} - }; - - }; - -} // namespace Opm - -#include "GasLiftRuntime_impl.hpp" - -#endif // OPM_GASLIFT_RUNTIME_HEADER_INCLUDED diff --git a/opm/simulators/wells/GasLiftRuntime_impl.hpp b/opm/simulators/wells/GasLiftRuntime_impl.hpp deleted file mode 100644 index 165e2a40f..000000000 --- a/opm/simulators/wells/GasLiftRuntime_impl.hpp +++ /dev/null @@ -1,816 +0,0 @@ -/* - Copyright 2020 Equinor ASA. - - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . -*/ - -#include -#include -#include -#include -#include -#include -#include - -#include -#include - -template -Opm::GasLiftRuntime:: -GasLiftRuntime( - const StdWell &std_well, - const Simulator &ebos_simulator, - const SummaryState &summary_state, - DeferredLogger &deferred_logger, - WellState &well_state, - const Well::ProductionControls &controls -) : - controls_{controls}, - deferred_logger_{deferred_logger}, - ebos_simulator_{ebos_simulator}, - potentials_(well_state.numPhases(),0.0), - std_well_{std_well}, - summary_state_{summary_state}, - well_state_{well_state}, - debug{false} // extra debugging output -{ - int well_index = this->std_well_.indexOfWell(); - const Well::ProducerCMode& control_mode - = well_state_.currentProductionControls()[well_index]; - if (control_mode != Well::ProducerCMode::THP) - throw std::logic_error("Bug in flow - invalid control mode detected\n"); - const Opm::Schedule& schedule = this->ebos_simulator_.vanguard().schedule(); - const int report_step_idx = this->ebos_simulator_.episodeIndex(); - auto ecl_well = this->std_well_.wellEcl(); - this->well_name_ = ecl_well.name(); - const GasLiftOpt& glo = schedule.glo(report_step_idx); - // NOTE: According to LIFTOPT, item 1: - // "Increment size for lift gas injection rate. Lift gas is - // allocated to individual wells in whole numbers of the increment - // size. If gas lift optimization is no longer required, it can be - // turned off by entering a zero or negative number." - // NOTE: This condition was checked in doGasLiftOptimize() in StandardWell - // so it can be assumed that increment_ > 0 - this->increment_ = glo.gaslift_increment(); - assert( this->increment_ > 0); - // NOTE: The manual (see LIFTOPT, item 2) does not mention - // any default value or restrictions on the economic gradient. - // TODO: The value of the gradient would most likely be a positive - // number. Should we warn or fail on a negative value? - // A negative value for the economic gradient would mean that - // the oil production is decreasing with increased liftgas - // injection (which seems strange) - this->eco_grad_ = glo.min_eco_gradient(); - auto& gl_well = glo.well(this->well_name_); - - if(useFixedAlq_(gl_well)) { - updateWellStateAlqFixedValue_(gl_well); - this->optimize_ = false; // lift gas supply is fixed - } - else { - setAlqMaxRate_(gl_well); - this->optimize_ = true; - } - computeInitialWellRates_(); - - if(this->optimize_) { - setAlqMinRate_(gl_well); - // NOTE: According to item 4 in WLIFTOPT, this value does not - // have to be positive. - // TODO: Does it make sense to have a negative value? - this->alpha_w_ = gl_well.weight_factor(); - if (this->alpha_w_ <= 0 ) { - displayWarning_("Nonpositive value for alpha_w ignored"); - this->alpha_w_ = 1.0; - } - - // NOTE: According to item 6 in WLIFTOPT: - // "If this value is greater than zero, the incremental gas rate will influence - // the calculation of the incremental gradient and may be used - // to discourage the allocation of lift gas to wells which produce more gas." - // TODO: Does this mean that we should ignore this value if it - // is negative? - this->alpha_g_ = gl_well.inc_weight_factor(); - - const auto& pu = std_well_.phaseUsage(); - this->oil_pos_ = pu.phase_pos[Oil]; - this->gas_pos_ = pu.phase_pos[Gas]; - this->water_pos_ = pu.phase_pos[Water]; - this->new_alq_ = this->orig_alq_; - // TODO: adhoc value.. Should we keep max_iterations_ as a safety measure - // or does it not make sense to have it? - this->max_iterations_ = 1000; - } -} - -/**************************************** - * Methods in alphabetical order - ****************************************/ - -template -void -Opm::GasLiftRuntime:: -computeInitialWellRates_() -{ - // get the alq value used for this well for the previous time step, or - // if gas lift optimization has not been applied to this well yet, the - // default value is used. - this->orig_alq_ = this->well_state_.getALQ(this->well_name_); - // NOTE: compute initial rates with current ALQ - this->std_well_.computeWellRatesWithThpAlqProd( - this->ebos_simulator_, this->summary_state_, this->deferred_logger_, - this->potentials_, this->orig_alq_); - -} - -template -void -Opm::GasLiftRuntime:: -computeWellRates_(double bhp, std::vector &potentials) -{ - this->std_well_.computeWellRatesWithBhp( - this->ebos_simulator_, bhp, potentials, this->deferred_logger_); - -} - -template -void -Opm::GasLiftRuntime:: -debugShowIterationInfo_(OptimizeState &state, double alq) -{ - const std::string msg = fmt::format("iteration {}, ALQ = {}", state.it, alq); - this->displayDebugMessage_(msg); -} - -template -void -Opm::GasLiftRuntime:: -debugShowStartIteration_(double alq, bool increase) -{ - const std::string msg = - fmt::format("starting {} iteration, ALQ = {}, oilrate = {}", - (increase ? "increase" : "decrease"), - alq, - -this->potentials_[this->oil_pos_]); - this->displayDebugMessage_(msg); -} - -template -void -Opm::GasLiftRuntime:: -displayDebugMessage_(const std::string &msg) -{ - - const std::string message = fmt::format( - " GLIFT (DEBUG) : Well {} : {}", this->well_name_, msg); - this->deferred_logger_.info(message); -} - -template -void -Opm::GasLiftRuntime:: -displayWarning_(std::string msg) -{ - const std::string message = fmt::format( - "GAS LIFT OPTIMIZATION, WELL {} : {}", this->well_name_, msg); - this->deferred_logger_.warning("WARNING", message); -} - -// TODO: what if the gas_rate_target_ has been defaulted -// (i.e. value == 0, meaning: "No limit") but the -// oil_rate_target_ has not been defaulted ? -// If the new_oil_rate exceeds the oil_rate_target_ it is cut back, -// but the same cut-back will not happen for the new_gas_rate -// Seems like an inconsistency, since alq should in this -// case also be adjusted (to the smaller value that would -// give oil target rate) but then the gas rate would also be smaller? -// The effect of not reducing the gas rate (if it should be -// reduced?) is that a too large value is used in the -// computation of the economic gradient making the gradient -// smaller than it should be since the term appears in the denominator. -template -bool -Opm::GasLiftRuntime:: -getGasRateWithLimit_(double& new_rate, const std::vector &potentials) -{ - new_rate = -potentials[this->gas_pos_]; - bool limit = false; - if (this->controls_.hasControl(Well::ProducerCMode::GRAT)) { - auto target = this->controls_.gas_rate; - if (new_rate > target) { - new_rate = target; - limit = true; - } - } - return limit; -} - - -// NOTE: If the computed oil rate is larger than the target -// rate of the well, we reduce it to the target rate. This -// will make the economic gradient smaller than it would be -// if we did not reduce the rate, and it is less -// likely that the current gas lift increment will be -// accepted. -// TODO: If it still is accepted, we should ideally reduce the alq -// also since we also reduced the rate. This might involve -// some sort of iteration though.. -template -bool -Opm::GasLiftRuntime:: -getOilRateWithLimit_(double& new_rate, const std::vector &potentials) -{ - new_rate = -potentials[this->oil_pos_]; - bool limit = false; - if (this->controls_.hasControl(Well::ProducerCMode::ORAT)) { - auto target = this->controls_.oil_rate; - if (new_rate > target) { - new_rate = target; - limit = true; - } - } - if (this->controls_.hasControl(Well::ProducerCMode::LRAT)) { - auto target = this->controls_.liquid_rate; - auto oil_rate = -potentials[this->oil_pos_]; - auto water_rate = -potentials[this->water_pos_]; - auto liq_rate = oil_rate + water_rate; - if (liq_rate > target) { - new_rate = oil_rate * (target/liq_rate); - limit = true; - } - } - return limit; -} - -template -void -Opm::GasLiftRuntime:: -logSuccess_() -{ - - const std::string message = fmt::format( - "GLIFT, WELL {} {} ALQ from {} to {}", - this->well_name_, - ((this->new_alq_ > this->orig_alq_) ? "increased" : "decreased"), - this->orig_alq_, - this->new_alq_); - this->deferred_logger_.info(message); -} - -/* - At this point we know that this is a production well, and that its current - * control mode is THP. - * - * - We would like to check if it is possible to - * 1) increase the oil production by adding lift gas injection to the - * well, or if that is not possible, if we 2) should reduce the amount - * of lift gas injected due to a too small gain in oil production - * (with the current lift gas injection rate) - * - For 1) above, we should not add lift gas if it would cause an oil - * rate target to be exceeded, and for 2) we should not reduce the - * amount of liftgas injected below the minimum lift gas injection - * rate. - * - * NOTE: If reducing or adding lift-gas further would cause - * one of the well targets like ORAT, WRAT, GRAT, LRAT, CRAT, RESV, BHP, - * to become violated we should stop the lift gas optimization - * loop.. and then updateWellControls() will later (hopefully) switch the well's - * control mode from THP to the mode of the violated target. - * - * - Lift gas is added if it is economical viable depending on - * the ratio of oil gained compared to the amount of liftgas added. - * - * - Lift gas supply may be limited. - * - * - The current value of liftgas for the well is stored in the WellState object. - * - * - It is assumed that the oil production rate is concave function F - * of the amount of lift gas, such that it increases initially due to the - * reduced density of the mixture in the tubing. However, as the - * lift gas supply is increased further, friction pressure losses in the - * tubing become more important, and the production rate peaks and - * then starts to decrease. - * Since lift gas injection has a price, e.g. compression costs can - * be expressed as a cost per unit rate of lift gas injection, - * it must be balanced against the value of the extra amount of - * oil produced. Thus there is a "minimum economic gradient" of oil - * production rate versus lift gas injection rate, at which the - * value of the extra amount of oil produced by a small increase in - * the lift gas injection rate is equal to the cost of supplying the - * extra amount of lift gas. The optimum lift gas injection rate is then somewhat - * lower than the peak value. - * - * Based on this assumption, we know that if the gradient (derivative) of F is - * positive, but greater than the economic gradient (assuming the - * economic gradient is positive), we should add - * lift gas. On the other hand, if the gradient of F is negative or - * if it is positive but smaller than the economic gradient, the amount - * of lift gas injected should be decreased. - * - */ -template -void -Opm::GasLiftRuntime:: -runOptimize() -{ - if (this->optimize_) { - if (!tryIncreaseLiftGas_()) { - if (!tryDecreaseLiftGas_()) { - return; - } - } - logSuccess_(); - this->well_state_.setALQ(this->well_name_, this->new_alq_); - } - // NOTE: In addition to the new ALQ value, we also implicitly - // return this->potentials_ -} - -// INPUT: -// - increase (boolean) : -// - true : try increase the lift gas supply, -// - false : try decrease lift gas supply. -// -// OUTPUT: -// - return value: success (true/false) -// - potentials_ : updated well potentials if success -// - new_alq_ : updated alq if success -template -bool -Opm::GasLiftRuntime:: -runOptimizeLoop_(bool increase) -{ - auto cur_potentials = this->potentials_; // make copy, since we may fail.. - auto oil_rate = -cur_potentials[this->oil_pos_]; - auto gas_rate = -cur_potentials[this->gas_pos_]; - bool success = false; // did we succeed to increase alq? - auto cur_alq = this->orig_alq_; - if (cur_alq == 0 && !increase) // we don't decrease alq below zero - return false; - - auto alq = cur_alq; - OptimizeState state {*this, increase}; - if (this->debug) debugShowStartIteration_(alq, increase); - while (!state.stop_iteration && (++state.it <= this->max_iterations_)) { - if (state.checkWellRatesViolated(cur_potentials)) break; - if (state.checkAlqOutsideLimits(alq, oil_rate)) break; - alq = state.addOrSubtractAlqIncrement(alq); - if(this->debug) debugShowIterationInfo_(state, alq); - if (!state.computeBhpAtThpLimit(alq)) break; - // NOTE: if BHP is below limit, we set state.stop_iteration = true - auto bhp = state.getBhpWithLimit(); - computeWellRates_(bhp, cur_potentials); - double new_oil_rate = 0.0; - bool oil_is_limited = getOilRateWithLimit_(new_oil_rate, cur_potentials); - if (oil_is_limited && !increase) //if oil is limited we do not want to decrease - break; - double new_gas_rate = 0.0; - bool gas_is_limited = getGasRateWithLimit_(new_gas_rate, cur_potentials); - if (gas_is_limited && increase) // if gas is limited we do not want to increase - break; - auto gradient = state.calcGradient( - oil_rate, new_oil_rate, gas_rate, new_gas_rate); - if (state.checkEcoGradient(gradient)) { - if (state.it == 1) { - break; - } - else { - state.stop_iteration = true; - } - } - cur_alq = alq; - success = true; - oil_rate = new_oil_rate; - gas_rate = new_gas_rate; - } - if (state.it > this->max_iterations_) { - warnMaxIterationsExceeded_(); - } - if (success) { - this->potentials_ = cur_potentials; - this->new_alq_ = cur_alq; - } - return success; -} - -template -void -Opm::GasLiftRuntime:: -setAlqMaxRate_(const GasLiftOpt::Well &well) -{ - auto& max_alq_optional = well.max_rate(); - if (max_alq_optional) { - // NOTE: To prevent extrapolation of the VFP tables, any value - // entered here must not exceed the largest ALQ value in the well's VFP table. - this->max_alq_ = *max_alq_optional; - } - else { // i.e. WLIFTOPT, item 3 has been defaulted - // According to the Eclipse manual for WLIFTOPT, item 3: - // The default value should be set to the largest ALQ - // value in the well's VFP table - const auto& table = std_well_.vfp_properties_->getProd()->getTable( - this->controls_.vfp_table_number); - const auto& alq_values = table.getALQAxis(); - // Assume the alq_values are sorted in ascending order, so - // the last item should be the largest value: - this->max_alq_ = alq_values.back(); - } -} - -template -void -Opm::GasLiftRuntime:: -setAlqMinRate_(const GasLiftOpt::Well &well) -{ - // NOTE: According to WLIFTOPT item 5 : - // if min_rate() is negative, it means: allocate at least enough lift gas - // to enable the well to flow - // NOTE: "to enable the well to flow" : How to interpret this? - // We choose to interpret it to mean a positive oil rate as returned from - // - // computeWellRates_(bhp, cur_potentials); - // - // So even if the well is producing gas, if the oil rate is zero - // we say that the "well is not flowing". - // - // Note that if WECON item 2 is set, the well can be shut off - // before the flow rate reaches zero. Also, - // if bhp drops below the bhp lower limit, the well might switch to bhp - // control before the oil rate becomes zero. - - this->min_alq_ = well.min_rate(); - if (this->min_alq_ > 0) { - if (this->min_alq_ >= this->max_alq_) { - // NOTE: We reset the value to a negative value. - // negative value means: Allocate at least enough lift gas - // to allow the well to flow. - // TODO: Consider other options for resetting the value.. - this->min_alq_ = -1; - displayWarning_("Minimum ALQ value is larger than maximum ALQ value!" - " Resetting value."); - } - } -} - -template -bool -Opm::GasLiftRuntime:: -tryDecreaseLiftGas_() -{ - return runOptimizeLoop_(/*increase=*/ false); -} - -template -bool -Opm::GasLiftRuntime:: -tryIncreaseLiftGas_() -{ - return runOptimizeLoop_(/*increase=*/ true); -} - - -// Called when we should use a fixed ALQ value -template -void -Opm::GasLiftRuntime:: -updateWellStateAlqFixedValue_(const GasLiftOpt::Well &well) -{ - auto& max_alq_optional = well.max_rate(); - if (max_alq_optional) { - // According to WLIFTOPT, item 3: - // If item 2 is NO, then item 3 is regarded as the fixed - // lift gas injection rate for the well. - auto new_alq = *max_alq_optional; - this->well_state_.setALQ(this->well_name_, new_alq); - } - // else { - // // If item 3 is defaulted, the lift gas rate remains - // // unchanged at its current value. - //} - -} - -// Determine if we should use a fixed ALQ value. -// -// From the manual for WLIFTOPT, item 2: -// Is the well's lift gas injection rate to be calculated by the -// optimization facility? -// - YES : The well's lift gas injection rate is calculated by the -// optimization facility. -// - NO : The well's lift gas injection rate remains fixed at a -// value that can be set either in Item 3 of this keyword, or in -// Item 12 of keyword WCONPROD, or with keyword WELTARG. -template -bool -Opm::GasLiftRuntime:: -useFixedAlq_(const GasLiftOpt::Well &well) -{ - auto wliftopt_item2 = well.use_glo(); - if (wliftopt_item2) { - return false; - } - else { - // auto& max_alq_optional = well.max_rate(); - // if (max_alq_optional) { - // According to WLIFTOPT, item 3: - // If item 2 is NO, then item 3 is regarded as the fixed - // lift gas injection rate for the well. - // } - // else { - // If item 3 is defaulted, the lift gas rate remains - // unchanged at its current value. - // } - return true; - } -} - -template -void -Opm::GasLiftRuntime:: -warnMaxIterationsExceeded_() -{ - std::ostringstream ss; - ss << "Max iterations (" << this->max_iterations_ << ") exceeded in " - << "gas lift optimization for well " << this->well_name_; - deferred_logger_.warning("MAX_ITERATIONS_EXCEEDED", ss.str()); -} - -/**************************************** - * Methods declared in OptimizeState - ****************************************/ - -template -double -Opm::GasLiftRuntime::OptimizeState:: -addOrSubtractAlqIncrement(double alq) -{ - if (this->increase) { - alq += this->parent.increment_; - // NOTE: if max_alq_ was defaulted in WCONPROD, item 3, it has - // already been set to the largest value in the VFP table in - // the contructor of GasLiftRuntime - if (alq > this->parent.max_alq_) alq = this->parent.max_alq_; - } - else { - alq -= this->parent.increment_; - if (this->parent.min_alq_ > 0) { - if (alq < this->parent.min_alq_) alq = this->parent.min_alq_; - } - } - return alq; -} - -template -double -Opm::GasLiftRuntime::OptimizeState:: -calcGradient(double oil_rate, double new_oil_rate, double gas_rate, double new_gas_rate) -{ - auto dqo = new_oil_rate - oil_rate; - auto dqg = new_gas_rate - gas_rate; - // TODO: Should we do any error checks on the calculation of the - // gradient? - // NOTE: The eclipse techincal description (chapter 25) says: - // "The gas rate term in the denominator is subject to the - // constraint alpha_g_ * dqg >= 0.0" - auto gradient = (this->parent.alpha_w_ * dqo) / - (this->parent.increment_ + this->parent.alpha_g_*dqg); - return gradient; -} - -// NOTE: According to WLIFTOPT item 5 : -// if min_rate() is negative, it means: allocate at least enough lift gas -// to enable the well to flow -// We will interpret this as (see discussion above GasLiftRuntime() -// in this file): Allocate at least the amount of lift gas needed to -// get a positive oil production rate. -template -bool -Opm::GasLiftRuntime::OptimizeState:: -checkAlqOutsideLimits(double alq, double oil_rate) -{ - std::ostringstream ss; - bool result = false; - - if (this->increase) { - if (alq >= this->parent.max_alq_) { - ss << "ALQ >= " << this->parent.max_alq_ << " (max limit), " - << "stopping iteration"; - result = true; - } - else { - // NOTE: A negative min_alq_ means: allocate at least enough lift gas - // to enable the well to flow, see WCONPROD item 5. - if (this->parent.min_alq_ < 0) { - result = false; - } - else { - // NOTE: checking for a lower limit should not be necessary - // when increasing alq.. so this is just to catch an - // illegal state at an early point. - if (alq < this->parent.min_alq_ ) { - warn_("unexpected: alq below lower limit when trying to " - "increase lift gas. aborting iteration."); - result = true; - } - else { - result = false; - } - } - } - } - else { // we are decreasing lift gas - // NOTE: A negative min_alq_ means: allocate at least enough lift gas - // to enable the well to flow, see WCONPROD item 5. - if (this->parent.min_alq_ < 0) { - // If the oil rate is already zero or negative (non-flowing well) - // we assume we will not be able to increase it by decreasing the lift gas - if ( oil_rate <= 0 ) { - ss << "Oil rate ( " << oil_rate << " ) <= 0 when decreasing lift gas. " - << "We will not be able to make this well flowing by decreasing " - << "lift gas, stopping iteration."; - result = true; - } - else { - result = false; - } - } - else { - if (alq <= this->parent.min_alq_ ) { - ss << "ALQ <= " << this->parent.min_alq_ << " (min limit), " - << "stopping iteration"; - result = true; - } - else { - // NOTE: checking for an upper limit should not be necessary - // when decreasing alq.. so this is just to catch an - // illegal state at an early point. - if (alq >= this->parent.max_alq_) { - warn_( "unexpected: alq above upper limit when trying to " - "decrease lift gas. aborting iteration."); - result = true; - } - else { - result = false; - } - } - } - } - if (this->parent.debug) this->parent.displayDebugMessage_(ss.str()); - return result; -} - -template -bool -Opm::GasLiftRuntime::OptimizeState:: -checkEcoGradient(double gradient) -{ - std::ostringstream ss; - bool result = false; - - if (this->parent.debug) { - ss << "checking gradient: " << gradient; - } - if (this->increase) { - if (this->parent.debug) ss << " <= " << this->parent.eco_grad_ << " --> "; - if (gradient <= this->parent.eco_grad_) { - if (this->parent.debug) ss << "true"; - result = true; - } - else { - if (this->parent.debug) ss << "false"; - } - } - else { // decreasing lift gas - if (this->parent.debug) ss << " >= " << this->parent.eco_grad_ << " --> "; - if (gradient >= this->parent.eco_grad_) { - if (this->parent.debug) ss << "true"; - result = true; - } - else { - if (this->parent.debug) ss << "false"; - } - } - if (this->parent.debug) this->parent.displayDebugMessage_(ss.str()); - return result; -} - -template -bool -Opm::GasLiftRuntime::OptimizeState:: -checkRate(double rate, double limit, const std::string rate_str) -{ - - if (limit < rate ) { - if (this->parent.debug) { - const std::string msg = fmt::format( - "iteration {} : rate {} exceeds target rate {}. Stopping iteration", - this->it, rate_str, rate, limit); - this->parent.displayDebugMessage_(msg); - } - return true; - } - return false; -} - -template -bool -Opm::GasLiftRuntime::OptimizeState:: -checkWellRatesViolated(std::vector &potentials) -{ - if (this->parent.controls_.hasControl(Well::ProducerCMode::ORAT)) { - auto oil_rate = -potentials[this->parent.oil_pos_]; - if (this->checkRate(oil_rate, this->parent.controls_.oil_rate, "oil")) - return true; - } - if (this->parent.controls_.hasControl(Well::ProducerCMode::WRAT)) { - auto water_rate = -potentials[this->parent.water_pos_]; - if (this->checkRate(water_rate, this->parent.controls_.water_rate, "water")) - return true; - } - if (this->parent.controls_.hasControl(Well::ProducerCMode::GRAT)) { - auto gas_rate = -potentials[this->parent.gas_pos_]; - if (this->checkRate(gas_rate, this->parent.controls_.gas_rate, "gas")) - return true; - } - if (this->parent.controls_.hasControl(Well::ProducerCMode::LRAT)) { - auto oil_rate = -potentials[this->parent.oil_pos_]; - auto water_rate = -potentials[this->parent.water_pos_]; - auto liq_rate = oil_rate + water_rate; - if (this->checkRate(liq_rate, this->parent.controls_.liquid_rate, "liquid")) - return true; - } - // TODO: Also check RESV, see checkIndividualContraints() in - // WellInterface_impl.hpp - // TODO: Check group contraints? - - return false; -} - - - -template -bool -Opm::GasLiftRuntime::OptimizeState:: -computeBhpAtThpLimit(double alq) -{ - auto bhp_at_thp_limit = this->parent.std_well_.computeBhpAtThpLimitProdWithAlq( - this->parent.ebos_simulator_, - this->parent.summary_state_, - this->parent.deferred_logger_, - alq); - if (!bhp_at_thp_limit) { - const std::string msg = fmt::format( - "Failed in getting converged bhp potential for well {}", - this->parent.well_name_); - this->parent.deferred_logger_.warning( - "FAILURE_GETTING_CONVERGED_POTENTIAL", msg); - return false; - } - this->bhp = *bhp_at_thp_limit; - return true; -} - -// NOTE: When calculating the gradient, determine what the well would produce if -// the lift gas injection rate were increased by one increment. The -// production rates are adjusted if necessary to obey -// any rate or BHP limits that the well may be subject to. From this -// information, calculate the well's "weighted incremental -// gradient" -// -// TODO: What does it mean to "adjust the production rates" given a -// BHP limit? -// -template -double -Opm::GasLiftRuntime::OptimizeState:: -getBhpWithLimit() -{ - auto bhp_update = this->bhp; - if (this->parent.controls_.hasControl(Well::ProducerCMode::BHP)) { - auto limit = this->parent.controls_.bhp_limit; - // TODO: is it possible that bhp falls below the limit when - // adding lift gas? I.e. if this->increase == true.. - if (this->bhp < limit) { - // TODO: we keep the current alq, but it should probably - // be adjusted since we changed computed bhp. But how? - bhp_update = limit; - // Stop iteration, but first check the economic gradient - // with the bhp_update. If the gradient looks OK (see the - // main optimize loop) we keep the current ALQ value. - this->stop_iteration = true; - } - } - return bhp_update; -} diff --git a/opm/simulators/wells/GasLiftSingleWell.hpp b/opm/simulators/wells/GasLiftSingleWell.hpp new file mode 100644 index 000000000..1526a0b36 --- /dev/null +++ b/opm/simulators/wells/GasLiftSingleWell.hpp @@ -0,0 +1,235 @@ +/* + Copyright 2020 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_GASLIFT_SINGLE_WELL_HEADER_INCLUDED +#define OPM_GASLIFT_SINGLE_WELL_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +// NOTE: StandardWell.hpp includes ourself (GasLiftSingleWell.hpp), so we need +// to forward declare StandardWell for it to be defined in this file. +namespace Opm { + template class StandardWell; +} +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + template + class GasLiftSingleWell + { + using Simulator = GetPropType; + using WellState = WellStateFullyImplicitBlackoil; + using StdWell = Opm::StandardWell; + using GLiftWellState = Opm::GasLiftWellState; + // TODO: same definition with WellInterface, and + // WellStateFullyImplicitBlackoil eventually they should go + // to a common header file. + static const int Water = BlackoilPhases::Aqua; + static const int Oil = BlackoilPhases::Liquid; + static const int Gas = BlackoilPhases::Vapour; + static constexpr double ALQ_EPSILON = 1e-8; + struct OptimizeState; + class Stage2State; + public: + GasLiftSingleWell( + const StdWell &std_well, + const Simulator &ebos_simulator, + const SummaryState &summary_state, + DeferredLogger &deferred_logger, + WellState &well_state + ); + struct GradInfo; + std::optional calcIncOrDecGradient( + double oil_rate, double gas_rate, double alq, bool increase) const; + std::pair getStage2Rates(); + const WellInterface &getStdWell() const { return std_well_; } + bool hasStage2Rates(); + std::unique_ptr runOptimize(); + const std::string& name() const {return well_name_; } + void updateStage2State(const GradInfo &gi, bool increase); + + struct GradInfo + { + GradInfo() { } + + GradInfo(double grad_, double new_oil_rate_, bool oil_is_limited_, + double new_gas_rate_, bool gas_is_limited_, + double alq_, bool alq_is_limited_) : + grad{grad_}, + new_oil_rate{new_oil_rate_}, + oil_is_limited{oil_is_limited_}, + new_gas_rate{new_gas_rate_}, + gas_is_limited{gas_is_limited_}, + alq{alq_}, + alq_is_limited{alq_is_limited_} {} + double grad; + double new_oil_rate; + bool oil_is_limited; + double new_gas_rate; + bool gas_is_limited; + double alq; + bool alq_is_limited; + }; + + private: + std::pair, bool> addOrSubtractAlqIncrement_( + 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 checkInitialALQmodified_(double alq, double initial_alq) const; + bool checkWellRatesViolated_( + std::vector &potentials, + const std::function &callback, + bool increase + ); + std::optional computeBhpAtThpLimit_(double alq) const; + bool computeInitialWellRates_(std::vector &potentials); + void computeWellRates_( + double bhp, std::vector &potentials, bool debug_output=true) const; + 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 debugShowStartIteration_(double alq, bool increase, double oil_rate); + void debugShowTargets_(); + void displayDebugMessage_(const std::string &msg) const; + void displayWarning_(std::string warning); + std::pair getBhpWithLimit_(double bhp) const; + std::pair getGasRateWithLimit_( + const std::vector &potentials) const; + std::tuple getInitialRatesWithLimit_( + const std::vector &potentials); + std::pair getOilRateWithLimit_( + const std::vector &potentials) const; + std::tuple + increaseALQtoPositiveOilRate_(double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials); + std::tuple + increaseALQtoMinALQ_(double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials); + void logSuccess_(double alq); + std::tuple + reduceALQtoOilTarget_(double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials); + + std::unique_ptr runOptimize1_(); + std::unique_ptr runOptimize2_(); + std::unique_ptr runOptimizeLoop_(bool increase); + void setAlqMaxRate_(const GasLiftOpt::Well &well); + void setAlqMinRate_(const GasLiftOpt::Well &well); + std::unique_ptr tryIncreaseLiftGas_(); + std::unique_ptr tryDecreaseLiftGas_(); + void updateWellStateAlqFixedValue_(const GasLiftOpt::Well &well); + bool useFixedAlq_(const GasLiftOpt::Well &well); + void warnMaxIterationsExceeded_(); + + DeferredLogger &deferred_logger_; + const Simulator &ebos_simulator_; + const StdWell &std_well_; + const SummaryState &summary_state_; + WellState &well_state_; + std::string well_name_; + const Well &ecl_well_; + const Well::ProductionControls controls_; + int num_phases_; + bool debug; // extra debug output + bool debug_limit_increase_decrease_; + bool debug_abort_if_decrease_and_oil_is_limited_ = false; + bool debug_abort_if_increase_and_gas_is_limited_ = false; + + double alpha_w_; + double alpha_g_; + double eco_grad_; + int gas_pos_; + bool has_run_init_ = false; + double increment_; + double max_alq_; + int max_iterations_; + double min_alq_; + int oil_pos_; + bool optimize_; + double orig_alq_; + int water_pos_; + + struct OptimizeState + { + OptimizeState( GasLiftSingleWell &parent_, bool increase_ ) : + parent{parent_}, + increase{increase_}, + it{0}, + stop_iteration{false}, + bhp{-1} + {} + + GasLiftSingleWell &parent; + bool increase; + int it; + bool stop_iteration; + double bhp; + + std::pair,bool> addOrSubtractAlqIncrement(double alq); + double calcEcoGradient(double oil_rate, double new_oil_rate, + double gas_rate, double new_gas_rate); + bool checkAlqOutsideLimits(double alq, double oil_rate); + bool checkEcoGradient(double gradient); + bool checkNegativeOilRate(double oil_rate); + bool checkOilRateExceedsTarget(double oil_rate); + bool checkRate(double rate, double limit, const std::string &rate_str) const; + bool checkWellRatesViolated(std::vector &potentials); + bool computeBhpAtThpLimit(double alq); + void debugShowIterationInfo(double alq); + double getBhpWithLimit(); + void warn_(std::string msg) {parent.displayWarning_(msg);} + }; + + }; + +#include "GasLiftSingleWell_impl.hpp" + +} // namespace Opm + + +#endif // OPM_GASLIFT_SINGLE_WELL_HEADER_INCLUDED diff --git a/opm/simulators/wells/GasLiftSingleWell_impl.hpp b/opm/simulators/wells/GasLiftSingleWell_impl.hpp new file mode 100644 index 000000000..929932e0c --- /dev/null +++ b/opm/simulators/wells/GasLiftSingleWell_impl.hpp @@ -0,0 +1,1370 @@ +/* + Copyright 2020 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +template +GasLiftSingleWell:: +GasLiftSingleWell( + const StdWell &std_well, + const Simulator &ebos_simulator, + const SummaryState &summary_state, + DeferredLogger &deferred_logger, + WellState &well_state +) : + deferred_logger_{deferred_logger}, + ebos_simulator_{ebos_simulator}, + std_well_{std_well}, + summary_state_{summary_state}, + well_state_{well_state}, + ecl_well_{std_well_.wellEcl()}, + controls_{ecl_well_.productionControls(summary_state_)}, + num_phases_{well_state_.numPhases()}, + debug{false}, // extra debugging output + debug_limit_increase_decrease_{false} +{ + const Schedule& schedule = this->ebos_simulator_.vanguard().schedule(); + const int report_step_idx = this->ebos_simulator_.episodeIndex(); + this->well_name_ = ecl_well_.name(); + const GasLiftOpt& glo = schedule.glo(report_step_idx); + // NOTE: According to LIFTOPT, item 1: + // "Increment size for lift gas injection rate. Lift gas is + // allocated to individual wells in whole numbers of the increment + // size. If gas lift optimization is no longer required, it can be + // turned off by entering a zero or negative number." + // NOTE: This condition was checked in doGasLiftOptimize() in StandardWell + // so it can be assumed that increment_ > 0 + this->increment_ = glo.gaslift_increment(); + assert( this->increment_ > 0); + // NOTE: The manual (see LIFTOPT, item 2) does not mention + // any default value or restrictions on the economic gradient. + // TODO: The value of the gradient would most likely be a positive + // number. Should we warn or fail on a negative value? + // A negative value for the economic gradient would mean that + // the oil production is decreasing with increased liftgas + // injection (which seems strange) + this->eco_grad_ = glo.min_eco_gradient(); + auto& gl_well = glo.well(this->well_name_); + + if(useFixedAlq_(gl_well)) { + updateWellStateAlqFixedValue_(gl_well); + this->optimize_ = false; // lift gas supply is fixed + } + else { + setAlqMaxRate_(gl_well); + this->optimize_ = true; + } + const auto& pu = std_well_.phaseUsage(); + this->oil_pos_ = pu.phase_pos[Oil]; + this->gas_pos_ = pu.phase_pos[Gas]; + this->water_pos_ = pu.phase_pos[Water]; + // get the alq value used for this well for the previous iteration (a + // nonlinear iteration in assemble() in BlackoilWellModel). + // If gas lift optimization has not been applied to this well yet, the + // default value is used. + this->orig_alq_ = this->well_state_.getALQ(this->well_name_); + if(this->optimize_) { + setAlqMinRate_(gl_well); + // NOTE: According to item 4 in WLIFTOPT, this value does not + // have to be positive. + // TODO: Does it make sense to have a negative value? + this->alpha_w_ = gl_well.weight_factor(); + if (this->alpha_w_ <= 0 ) { + displayWarning_("Nonpositive value for alpha_w ignored"); + this->alpha_w_ = 1.0; + } + + // NOTE: According to item 6 in WLIFTOPT: + // "If this value is greater than zero, the incremental gas rate will influence + // the calculation of the incremental gradient and may be used + // to discourage the allocation of lift gas to wells which produce more gas." + // TODO: Does this mean that we should ignore this value if it + // is negative? + this->alpha_g_ = gl_well.inc_weight_factor(); + + // TODO: adhoc value.. Should we keep max_iterations_ as a safety measure + // or does it not make sense to have it? + this->max_iterations_ = 1000; + } +} + +/**************************************** + * Public methods in alphabetical order + ****************************************/ + +// NOTE: Used from GasLiftStage2 +template +std::optional::GradInfo> +GasLiftSingleWell:: +calcIncOrDecGradient(double oil_rate, double gas_rate, double alq, bool increase) const +{ + auto [new_alq_opt, alq_is_limited] = addOrSubtractAlqIncrement_(alq, increase); + // TODO: What to do if ALQ is limited and new_alq != alq? + if (!new_alq_opt) + return std::nullopt; + double new_alq = *new_alq_opt; + if (auto bhp = computeBhpAtThpLimit_(new_alq)) { + auto [new_bhp, bhp_is_limited] = getBhpWithLimit_(*bhp); + // TODO: What to do if BHP is limited? + std::vector potentials(this->num_phases_, 0.0); + computeWellRates_(new_bhp, potentials); + auto [new_oil_rate, oil_is_limited] = getOilRateWithLimit_(potentials); + auto [new_gas_rate, gas_is_limited] = getGasRateWithLimit_(potentials); + if (!increase && new_oil_rate < 0 ) { + return std::nullopt; + } + auto grad = calcEcoGradient_( + oil_rate, new_oil_rate, gas_rate, new_gas_rate, increase); + return GradInfo(grad, new_oil_rate, oil_is_limited, + new_gas_rate, gas_is_limited, new_alq, alq_is_limited); + } + else { + return std::nullopt; + } +} + +/* - At this point we know that this is a production well, and that its current + * control mode is THP. + * + * - We would like to check if it is possible to + * 1) increase the oil production by adding lift gas injection to the + * well, or if that is not possible, if we 2) should reduce the amount + * of lift gas injected due to a too small gain in oil production + * (with the current lift gas injection rate) + * - For 1) above, we should not add lift gas if it would cause an oil + * rate target to be exceeded, and for 2) we should not reduce the + * amount of liftgas injected below the minimum lift gas injection + * rate. + * + * NOTE: If reducing or adding lift-gas further would cause + * one of the well targets like ORAT, WRAT, GRAT, LRAT, CRAT, RESV, BHP, + * to become violated we should stop the lift gas optimization + * loop.. and then updateWellControls() will later (hopefully) switch the well's + * control mode from THP to the mode of the violated target. + * + * - Lift gas is added if it is economical viable depending on + * the ratio of oil gained compared to the amount of liftgas added. + * + * - Lift gas supply may be limited. + * + * - The current value of liftgas for the well is stored in the WellState object. + * + * - It is assumed that the oil production rate is concave function F + * of the amount of lift gas, such that it increases initially due to the + * reduced density of the mixture in the tubing. However, as the + * lift gas supply is increased further, friction pressure losses in the + * tubing become more important, and the production rate peaks and + * then starts to decrease. + * Since lift gas injection has a price, e.g. compression costs can + * be expressed as a cost per unit rate of lift gas injection, + * it must be balanced against the value of the extra amount of + * oil produced. Thus there is a "minimum economic gradient" of oil + * production rate versus lift gas injection rate, at which the + * value of the extra amount of oil produced by a small increase in + * the lift gas injection rate is equal to the cost of supplying the + * extra amount of lift gas. The optimum lift gas injection rate is then somewhat + * lower than the peak value. + * + * Based on this assumption, we know that if the gradient (derivative) of F is + * positive, but greater than the economic gradient (assuming the + * economic gradient is positive), we should add + * lift gas. On the other hand, if the gradient of F is negative or + * if it is positive but smaller than the economic gradient, the amount + * of lift gas injected should be decreased. + * + */ +template +std::unique_ptr::GLiftWellState> +GasLiftSingleWell:: +runOptimize() +{ + std::unique_ptr state; + if (this->optimize_) { + if (this->debug_limit_increase_decrease_) { + state = runOptimize1_(); + } + else { + state = runOptimize2_(); + } + if (state) { + if (state->increase()) { + double alq = state->alq(); + if (this->debug) + logSuccess_(alq); + this->well_state_.setALQ(this->well_name_, alq); + } + } + } + return state; +} + +template +std::unique_ptr::GLiftWellState> +GasLiftSingleWell:: +runOptimize2_() +{ + std::unique_ptr state; + state = tryIncreaseLiftGas_(); + if (!state || !(state->alqChanged())) { + state = tryDecreaseLiftGas_(); + } + return state; +} + +template +std::unique_ptr::GLiftWellState> +GasLiftSingleWell:: +runOptimize1_() +{ + std::unique_ptr state; + auto inc_count = this->well_state_.gliftGetAlqIncreaseCount(this->well_name_); + auto dec_count = this->well_state_.gliftGetAlqDecreaseCount(this->well_name_); + if (dec_count == 0 && inc_count == 0) { + state = tryIncreaseLiftGas_(); + if (!state && !state->alqChanged()) { + state = tryDecreaseLiftGas_(); + } + } + else if (dec_count == 0) { + assert(inc_count > 0); + state = tryIncreaseLiftGas_(); + } + else if (inc_count == 0) { + assert(dec_count > 0); + state = tryDecreaseLiftGas_(); + } + return state; +} + + +/**************************************** + * Private methods in alphabetical order + ****************************************/ + +template +std::pair, bool> +GasLiftSingleWell:: +addOrSubtractAlqIncrement_(double alq, bool increase) const +{ + bool limited = false; + double orig_alq = alq; + if (increase) { + alq += this->increment_; + // NOTE: if max_alq_ was defaulted in WLIFTOPT, item 3, it has + // already been set to the largest value in the VFP table in + // the contructor of GasLiftSingleWell + if (alq > this->max_alq_) { + alq = this->max_alq_; + limited = true; + } + } + else { // we are decreasing ALQ + alq -= this->increment_; + if (this->min_alq_ > 0) { + // According to WLIFTOPT item 5: If a positive value is + // specified (min_alq_), the well is allocated at least that amount + // of lift gas, unless the well is unable to flow with + // that rate of lift gas injection, or unless the well can + // already meet one of its own rate limits before + // receiving its minimum lift gas rate. + if (alq < this->min_alq_) { + alq = this->min_alq_; + limited = true; + } + } + else { + if (alq < 0) { + alq = 0.0; + limited = true; + } + } + } + std::optional alq_opt {alq}; + // If we were not able to change ALQ (to within rounding error), we + // return std::nullopt + if (limited && checkALQequal_(orig_alq,alq)) + alq_opt = std::nullopt; + + return {alq_opt, limited}; +} + +template +double +GasLiftSingleWell:: +calcEcoGradient_( + double oil_rate, double new_oil_rate, double gas_rate, + double new_gas_rate, bool increase) const +{ + auto dqo = new_oil_rate - oil_rate; + auto dqg = new_gas_rate - gas_rate; + // TODO: Should the gas rate term in the denominator be subject to the constraint: + // + // alpha_g_ * dqg >= 0.0 + // + // ? + auto gradient = (this->alpha_w_ * dqo) / (this->increment_ + this->alpha_g_*dqg); + // TODO: Should we do any error checks on the calculation of the + // gradient? + + if (!increase) gradient = -gradient; + return gradient; +} + +template +bool +GasLiftSingleWell:: +checkALQequal_(double alq1, double alq2) const +{ + return std::fabs(alq1-alq2) < (this->increment_*ALQ_EPSILON); +} + +template +bool +GasLiftSingleWell:: +checkInitialALQmodified_(double alq, double initial_alq) const +{ + if (checkALQequal_(alq,initial_alq)) { + return false; + } + else { + const std::string msg = fmt::format("initial ALQ changed from {} " + "to {} before iteration starts..", initial_alq, alq); + displayDebugMessage_(msg); + return true; + } +} + + +template +bool +GasLiftSingleWell:: +checkWellRatesViolated_( + std::vector &potentials, + const std::function &callback, + bool increase +) +{ + if (!increase) { + auto oil_rate = -potentials[this->oil_pos_]; + if (oil_rate < 0) { + // The well is not flowing, and it will(?) not help to reduce lift + // gas further. Note that this assumes that the oil rates drops with + // decreasing lift gas. + displayDebugMessage_("Negative oil rate detected while descreasing " + "lift gas. Stopping iteration."); + return true; + } + } + // TODO: the below checks could probably be skipped if we are decreasing + // lift gas (provided we can assume that rates declines monotonically with + // decreasing lift gas). + if (this->controls_.hasControl(Well::ProducerCMode::ORAT)) { + auto oil_rate = -potentials[this->oil_pos_]; + if (callback(oil_rate, this->controls_.oil_rate, "oil")) + return true; + } + if (this->controls_.hasControl(Well::ProducerCMode::WRAT)) { + auto water_rate = -potentials[this->water_pos_]; + if (callback(water_rate, this->controls_.water_rate, "water")) + return true; + } + if (this->controls_.hasControl(Well::ProducerCMode::GRAT)) { + auto gas_rate = -potentials[this->gas_pos_]; + if (callback(gas_rate, this->controls_.gas_rate, "gas")) + return true; + } + if (this->controls_.hasControl(Well::ProducerCMode::LRAT)) { + auto oil_rate = -potentials[this->oil_pos_]; + auto water_rate = -potentials[this->water_pos_]; + auto liq_rate = oil_rate + water_rate; + if (callback(liq_rate, this->controls_.liquid_rate, "liquid")) + return true; + } + // TODO: Also check RESV, see checkIndividualContraints() in + // WellInterface_impl.hpp + // TODO: Check group contraints? + + return false; +} + + +template +std::optional +GasLiftSingleWell:: +computeBhpAtThpLimit_(double alq) const +{ + auto bhp_at_thp_limit = this->std_well_.computeBhpAtThpLimitProdWithAlq( + this->ebos_simulator_, + this->summary_state_, + this->deferred_logger_, + alq); + 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); + bhp_at_thp_limit = this->controls_.bhp_limit; + } + //bhp_at_thp_limit = std::max(*bhp_at_thp_limit, this->controls_.bhp_limit); + } + else { + const std::string msg = fmt::format( + "Failed in getting converged bhp potential from thp limit (ALQ = {})", alq); + displayDebugMessage_(msg); + } + return bhp_at_thp_limit; +} + +template +bool +GasLiftSingleWell:: +computeInitialWellRates_(std::vector &potentials) +{ + if (auto bhp = computeBhpAtThpLimit_(this->orig_alq_); bhp) { + { + const std::string msg = fmt::format( + "computed initial bhp {} given thp limit and given alq {}", + *bhp, this->orig_alq_); + displayDebugMessage_(msg); + } + computeWellRates_(*bhp, potentials); + { + const std::string msg = fmt::format( + "computed initial well potentials given bhp, " + "oil: {}, gas: {}, water: {}", + -potentials[this->oil_pos_], + -potentials[this->gas_pos_], + -potentials[this->water_pos_]); + displayDebugMessage_(msg); + } + return true; + } + else { + displayDebugMessage_("Aborting optimization."); + return false; + } +} + +template +void +GasLiftSingleWell:: +computeWellRates_( + double bhp, std::vector &potentials, bool debug_output) const +{ + // NOTE: If we do not clear the potentials here, it will accumulate + // 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) { + const std::string msg = fmt::format("computed well potentials given bhp {}, " + "oil: {}, gas: {}, water: {}", bhp, + -potentials[this->oil_pos_], -potentials[this->gas_pos_], + -potentials[this->water_pos_]); + displayDebugMessage_(msg); + } +} + +template +void +GasLiftSingleWell:: +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 +{ + { + const std::string msg = fmt::format("calculating gradient: " + "new_oil_rate = {}, oil_rate = {}, grad = {}", new_oil_rate, oil_rate, grad); + displayDebugMessage_(msg); + } + if (grad < 0 ) { + const std::string msg = fmt::format("negative {} gradient detected ({}) : " + "alq: {}, new_alq: {}, " + "oil_rate: {}, new_oil_rate: {}, gas_rate: {}, new_gas_rate: {}", + (increase ? "incremental" : "decremental"), + grad, alq, new_alq, oil_rate, new_oil_rate, gas_rate, new_gas_rate); + displayDebugMessage_(msg); + } +} + +template +void +GasLiftSingleWell:: +debugShowBhpAlqTable_() +{ + double alq = 0.0; + const std::string fmt_fmt1 {"{:^12s} {:^12s} {:^12s} {:^12s}"}; + const std::string fmt_fmt2 {"{:>12.5g} {:>12.5g} {:>12.5g} {:>12.5g}"}; + const std::string header = fmt::format(fmt_fmt1, "ALQ", "BHP", "oil", "gas"); + displayDebugMessage_(header); + while (alq <= (this->max_alq_+this->increment_)) { + auto bhp_at_thp_limit = computeBhpAtThpLimit_(alq); + if (!bhp_at_thp_limit) { + const std::string msg = fmt::format("Failed to get converged potentials " + "for ALQ = {}. Skipping.", alq ); + displayDebugMessage_(msg); + } + else { + std::vector potentials(this->num_phases_, 0.0); + computeWellRates_(*bhp_at_thp_limit, potentials, /*debug_out=*/false); + auto oil_rate = -potentials[this->oil_pos_]; + auto gas_rate = -potentials[this->gas_pos_]; + const std::string msg = fmt::format( + fmt_fmt2, alq, *bhp_at_thp_limit, oil_rate, gas_rate); + displayDebugMessage_(msg); + } + alq += this->increment_; + } +} + +template +void +GasLiftSingleWell:: +debugShowStartIteration_(double alq, bool increase, double oil_rate) +{ + const std::string msg = + fmt::format("starting {} iteration, ALQ = {}, oilrate = {}", + (increase ? "increase" : "decrease"), + alq, oil_rate); + displayDebugMessage_(msg); +} + +template +void +GasLiftSingleWell:: +debugShowTargets_() +{ + if (this->controls_.hasControl(Well::ProducerCMode::ORAT)) { + auto target = this->controls_.oil_rate; + const std::string msg = fmt::format("has ORAT control with target {}", target); + displayDebugMessage_(msg); + } + if (this->controls_.hasControl(Well::ProducerCMode::GRAT)) { + auto target = this->controls_.gas_rate; + const std::string msg = fmt::format("has GRAT control with target {}", target); + displayDebugMessage_(msg); + } + if (this->controls_.hasControl(Well::ProducerCMode::LRAT)) { + auto target = this->controls_.liquid_rate; + const std::string msg = fmt::format("has LRAT control with target {}", target); + displayDebugMessage_(msg); + } +} + +template +void +GasLiftSingleWell:: +debugShowAlqIncreaseDecreaseCounts_() +{ + auto inc_count = this->well_state_.gliftGetAlqIncreaseCount(this->well_name_); + auto dec_count = this->well_state_.gliftGetAlqDecreaseCount(this->well_name_); + const std::string msg = + fmt::format("ALQ increase/decrease count : {}/{}", inc_count, dec_count); + displayDebugMessage_(msg); + +} + +template +void +GasLiftSingleWell:: +displayDebugMessage_(const std::string &msg) const +{ + + if (this->debug) { + const std::string message = fmt::format( + " GLIFT (DEBUG) : Well {} : {}", this->well_name_, msg); + this->deferred_logger_.info(message); + } +} + +template +void +GasLiftSingleWell:: +displayWarning_(std::string msg) +{ + const std::string message = fmt::format( + "GAS LIFT OPTIMIZATION, WELL {} : {}", this->well_name_, msg); + this->deferred_logger_.warning("WARNING", message); +} + +template +std::pair +GasLiftSingleWell:: +getBhpWithLimit_(double bhp) const +{ + bool limited = false; + if (this->controls_.hasControl(Well::ProducerCMode::BHP)) { + auto limit = this->controls_.bhp_limit; + if (bhp < limit) { + bhp = limit; + limited = true; + } + } + return {bhp, limited}; +} + +// TODO: what if the gas_rate_target_ has been defaulted +// (i.e. value == 0, meaning: "No limit") but the +// oil_rate_target_ has not been defaulted ? +// If the new_oil_rate exceeds the oil_rate_target_ it is cut back, +// but the same cut-back will not happen for the new_gas_rate +// Seems like an inconsistency, since alq should in this +// case also be adjusted (to the smaller value that would +// give oil target rate) but then the gas rate would also be smaller? +// The effect of not reducing the gas rate (if it should be +// reduced?) is that a too large value is used in the +// computation of the economic gradient making the gradient +// smaller than it should be since the term appears in the denominator. +template +std::pair +GasLiftSingleWell:: +getGasRateWithLimit_(const std::vector &potentials) const +{ + double new_rate = -potentials[this->gas_pos_]; + bool limit = false; + if (this->controls_.hasControl(Well::ProducerCMode::GRAT)) { + auto target = this->controls_.gas_rate; + if (new_rate > target) { + new_rate = target; + limit = true; + } + } + return { new_rate, limit}; +} + + +// NOTE: If the computed oil rate is larger than the target +// rate of the well, we reduce it to the target rate. This +// will make the economic gradient smaller than it would be +// if we did not reduce the rate, and it is less +// likely that the current gas lift increment will be +// accepted. +// TODO: If it still is accepted, we should ideally reduce the alq +// also since we also reduced the rate. This might involve +// some sort of iteration though.. +template +std::pair +GasLiftSingleWell:: +getOilRateWithLimit_(const std::vector &potentials) const +{ + double new_rate = -potentials[this->oil_pos_]; + if (this->controls_.hasControl(Well::ProducerCMode::ORAT)) { + auto target = this->controls_.oil_rate; + if (new_rate > target) { + const std::string msg = fmt::format("limiting oil rate to target: " + "computed rate: {}, target: {}", new_rate, target); + displayDebugMessage_(msg); + new_rate = target; + return { new_rate, /*limit=*/true}; + } + } + // TODO: What to do if both oil and lrat are violated? + // Currently oil rate violation takes precedence if both are + // violated, and the new rate is set to the oil rate target. + if (this->controls_.hasControl(Well::ProducerCMode::LRAT)) { + auto target = this->controls_.liquid_rate; + auto oil_rate = -potentials[this->oil_pos_]; + auto water_rate = -potentials[this->water_pos_]; + auto liq_rate = oil_rate + water_rate; + if (liq_rate > target) { + new_rate = oil_rate * (target/liq_rate); + const std::string msg = fmt::format( + "limiting oil rate due to LRAT target {}: " + "computed rate: {}, target: {}", target, oil_rate, new_rate); + displayDebugMessage_(msg); + return { new_rate, /*limit=*/true}; + } + } + return { new_rate, /*limit=*/false}; +} + +template +std::tuple +GasLiftSingleWell:: +getInitialRatesWithLimit_(const std::vector &potentials) +{ + auto [oil_rate, oil_is_limited] = getOilRateWithLimit_(potentials); + auto [gas_rate, gas_is_limited] = getGasRateWithLimit_(potentials); + + if (oil_is_limited) { + const std::string msg = fmt::format( + "initial oil rate was limited to: {}", oil_rate); + displayDebugMessage_(msg); + } + if (gas_is_limited) { + const std::string msg = fmt::format( + "initial gas rate was limited to: {}", gas_rate); + displayDebugMessage_(msg); + } + return std::make_tuple(oil_rate, gas_rate, oil_is_limited, gas_is_limited); +} + +template +std::tuple +GasLiftSingleWell:: +increaseALQtoPositiveOilRate_(double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials) +{ + bool stop_iteration = false; + double temp_alq = alq; + while(!stop_iteration) { + temp_alq += this->increment_; + if (temp_alq > this->max_alq_) break; + auto bhp_opt = computeBhpAtThpLimit_(temp_alq); + if (!bhp_opt) break; + alq = temp_alq; + auto [bhp, bhp_is_limited] = getBhpWithLimit_(*bhp_opt); + computeWellRates_(bhp, potentials); + oil_rate = -potentials[this->oil_pos_]; + if (oil_rate > 0) break; + } + std::tie(oil_rate, oil_is_limited) = getOilRateWithLimit_(potentials); + std::tie(gas_rate, gas_is_limited) = getGasRateWithLimit_(potentials); + return std::make_tuple(oil_rate, gas_rate, oil_is_limited, gas_is_limited, alq); +} + +template +std::tuple +GasLiftSingleWell:: +increaseALQtoMinALQ_(double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials) +{ + auto min_alq = this->min_alq_; + assert(min_alq >= 0); + assert(alq < min_alq); + assert(min_alq < this->max_alq_); + bool stop_iteration = false; + double temp_alq = alq; + while(!stop_iteration) { + temp_alq += this->increment_; + if (temp_alq >= min_alq) break; + auto bhp_opt = computeBhpAtThpLimit_(temp_alq); + if (!bhp_opt) break; + alq = temp_alq; + auto [bhp, bhp_is_limited] = getBhpWithLimit_(*bhp_opt); + computeWellRates_(bhp, potentials); + std::tie(oil_rate, oil_is_limited) = getOilRateWithLimit_(potentials); + std::tie(gas_rate, gas_is_limited) = getGasRateWithLimit_(potentials); + if (oil_is_limited || gas_is_limited) break; + } + return std::make_tuple(oil_rate, gas_rate, oil_is_limited, gas_is_limited, alq); +} + +template +void +GasLiftSingleWell:: +logSuccess_(double alq) +{ + const int iteration_idx = + this->ebos_simulator_.model().newtonMethod().numIterations(); + const std::string message = fmt::format( + "GLIFT, IT={}, WELL {} : {} ALQ from {} to {}", + iteration_idx, + this->well_name_, + ((alq > this->orig_alq_) ? "increased" : "decreased"), + this->orig_alq_, alq); + this->deferred_logger_.info(message); +} + +template +std::tuple +GasLiftSingleWell:: +reduceALQtoOilTarget_(double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials) +{ + displayDebugMessage_("Reducing ALQ to meet oil target before iteration starts.."); + double orig_oil_rate = oil_rate; + double orig_alq = alq; + // NOTE: This method should only be called if oil_is_limited, and hence + // we know that it has oil rate control + assert(this->controls_.hasControl(Well::ProducerCMode::ORAT)); + auto target = this->controls_.oil_rate; + bool stop_iteration = false; + double temp_alq = alq; + while(!stop_iteration) { + temp_alq -= this->increment_; + if (temp_alq <= 0) break; + auto bhp_opt = computeBhpAtThpLimit_(temp_alq); + if (!bhp_opt) break; + alq = temp_alq; + auto [bhp, bhp_is_limited] = getBhpWithLimit_(*bhp_opt); + computeWellRates_(bhp, potentials); + oil_rate = -potentials[this->oil_pos_]; + if (oil_rate < target) { + break; + } + } + std::tie(oil_rate, oil_is_limited) = getOilRateWithLimit_(potentials); + std::tie(gas_rate, gas_is_limited) = getGasRateWithLimit_(potentials); + if (this->debug) { + const std::string msg = fmt::format("Reduced oil rate from {} to {} and ALQ ", + "from {} to {}", orig_oil_rate, oil_rate, orig_alq, alq); + displayDebugMessage_(msg); + } + return std::make_tuple(oil_rate, gas_rate, oil_is_limited, gas_is_limited, alq); +} + +// INPUT: +// - increase (boolean) : +// - true : try increase the lift gas supply, +// - false : try decrease lift gas supply. +// +// OUTPUT: +// +// - return value: a new GLiftWellState or nullptr +// +template +std::unique_ptr::GLiftWellState> +GasLiftSingleWell:: +runOptimizeLoop_(bool increase) +{ + std::vector potentials(this->num_phases_, 0.0); + std::unique_ptr ret_value; // nullptr initially + if (!computeInitialWellRates_(potentials)) return ret_value; + bool alq_is_limited, oil_is_limited, gas_is_limited; + double oil_rate, gas_rate; + std::tie(oil_rate, gas_rate, oil_is_limited, gas_is_limited) = + getInitialRatesWithLimit_(potentials); + if (this->debug) debugShowBhpAlqTable_(); + if (this->debug) debugShowAlqIncreaseDecreaseCounts_(); + if (this->debug) debugShowTargets_(); + bool success = false; // did we succeed to increase alq? + auto cur_alq = this->orig_alq_; + auto temp_alq = cur_alq; + if (!increase && oil_is_limited) { + // NOTE: Try to decrease ALQ down to a value where the oil target is + // not exceeded. + // NOTE: This may reduce ALQ below the minimum value set in WLIFTOPT + // item 5. However, this is OK since the rate target is met and there + // is no point in using a higher ALQ value then. + std::tie(oil_rate, gas_rate, oil_is_limited, gas_is_limited, temp_alq) = + reduceALQtoOilTarget_( + temp_alq, oil_rate, gas_rate, oil_is_limited, gas_is_limited, potentials); + } + else { + if (increase && oil_rate < 0) { + // NOTE: Try to increase ALQ up to a value where oil_rate is positive + std::tie(oil_rate, gas_rate, oil_is_limited, gas_is_limited, temp_alq) = + increaseALQtoPositiveOilRate_(temp_alq, oil_rate, + gas_rate, oil_is_limited, gas_is_limited, potentials); + } + if (increase && (this->min_alq_> 0) && (temp_alq < this->min_alq_)) { + // NOTE: Try to increase ALQ up to the minimum limit without checking + // the economic gradient.. + std::tie(oil_rate, gas_rate, oil_is_limited, gas_is_limited, temp_alq) = + increaseALQtoMinALQ_(temp_alq, oil_rate, gas_rate, + oil_is_limited, gas_is_limited, potentials); + } + } + if (checkInitialALQmodified_(temp_alq, cur_alq)) { + cur_alq = temp_alq; + success = true; + } + OptimizeState state {*this, increase}; + if (this->debug) debugShowStartIteration_(temp_alq, increase, oil_rate); + while (!state.stop_iteration && (++state.it <= this->max_iterations_)) { + if (!increase && state.checkNegativeOilRate(oil_rate)) break; + if (state.checkWellRatesViolated(potentials)) break; + if (state.checkAlqOutsideLimits(temp_alq, oil_rate)) break; + std::optional alq_opt; + std::tie(alq_opt, alq_is_limited) + = state.addOrSubtractAlqIncrement(temp_alq); + if (!alq_opt) break; + temp_alq = *alq_opt; + if (this->debug) state.debugShowIterationInfo(temp_alq); + if (!state.computeBhpAtThpLimit(temp_alq)) break; + // NOTE: if BHP is below limit, we set state.stop_iteration = true + auto bhp = state.getBhpWithLimit(); + computeWellRates_(bhp, potentials); + auto [new_oil_rate, new_oil_is_limited] = getOilRateWithLimit_(potentials); +/* if (this->debug_abort_if_decrease_and_oil_is_limited_) { + if (oil_is_limited && !increase) { + // if oil is limited we do not want to decrease + displayDebugMessage_( + "decreasing ALQ and oil is limited -> aborting iteration"); + break; + } + } +*/ + auto [new_gas_rate, new_gas_is_limited] = getGasRateWithLimit_(potentials); +/* if (this->debug_abort_if_increase_and_gas_is_limited_) { + if (gas_is_limited && increase) { + // if gas is limited we do not want to increase + displayDebugMessage_( + "increasing ALQ and gas is limited -> aborting iteration"); + break; + } + } +*/ + auto gradient = state.calcEcoGradient( + oil_rate, new_oil_rate, gas_rate, new_gas_rate); + if (this->debug) + debugCheckNegativeGradient_( + gradient, cur_alq, temp_alq, oil_rate, new_oil_rate, + gas_rate, new_gas_rate, increase); + if (state.checkEcoGradient(gradient)) break; + cur_alq = temp_alq; + success = true; + oil_rate = new_oil_rate; + gas_rate = new_gas_rate; + oil_is_limited = new_oil_is_limited; + gas_is_limited = new_gas_is_limited; + } + if (state.it > this->max_iterations_) { + warnMaxIterationsExceeded_(); + } + std::optional increase_opt; + if (success) { + this->well_state_.gliftUpdateAlqIncreaseCount(this->well_name_, increase); + increase_opt = increase; + } + else { + increase_opt = std::nullopt; + } + ret_value = std::make_unique(oil_rate, oil_is_limited, + gas_rate, gas_is_limited, cur_alq, alq_is_limited, increase_opt); + return ret_value; +} + +template +void +GasLiftSingleWell:: +setAlqMaxRate_(const GasLiftOpt::Well &well) +{ + auto& max_alq_optional = well.max_rate(); + if (max_alq_optional) { + // NOTE: To prevent extrapolation of the VFP tables, any value + // entered here must not exceed the largest ALQ value in the well's VFP table. + this->max_alq_ = *max_alq_optional; + } + else { // i.e. WLIFTOPT, item 3 has been defaulted + // According to the manual for WLIFTOPT, item 3: + // The default value should be set to the largest ALQ + // value in the well's VFP table + const auto& table = std_well_.vfp_properties_->getProd()->getTable( + this->controls_.vfp_table_number); + const auto& alq_values = table.getALQAxis(); + // Assume the alq_values are sorted in ascending order, so + // the last item should be the largest value: + this->max_alq_ = alq_values.back(); + } +} + +template +void +GasLiftSingleWell:: +setAlqMinRate_(const GasLiftOpt::Well &well) +{ + // NOTE: According to WLIFTOPT item 5 : + // if min_rate() is negative, it means: allocate at least enough lift gas + // to enable the well to flow + // NOTE: "to enable the well to flow" : How to interpret this? + // We choose to interpret it to mean a positive oil rate as returned from + // + // computeWellRates_(bhp, cur_potentials); + // + // So even if the well is producing gas, if the oil rate is zero + // we say that the "well is not flowing". + // + // Note that if WECON item 2 is set, the well can be shut off + // before the flow rate reaches zero. Also, + // if bhp drops below the bhp lower limit, the well might switch to bhp + // control before the oil rate becomes zero. + + this->min_alq_ = well.min_rate(); + if (this->min_alq_ > 0) { + if (this->min_alq_ >= this->max_alq_) { + // NOTE: We reset the value to a negative value. + // negative value means: Allocate at least enough lift gas + // to allow the well to flow. + // TODO: Consider other options for resetting the value.. + this->min_alq_ = -1; + displayWarning_("Minimum ALQ value is larger than maximum ALQ value!" + " Resetting value."); + } + } +} + +template +std::unique_ptr::GLiftWellState> +GasLiftSingleWell:: +tryDecreaseLiftGas_() +{ + return runOptimizeLoop_(/*increase=*/ false); +} + +template +std::unique_ptr::GLiftWellState> +GasLiftSingleWell:: +tryIncreaseLiftGas_() +{ + return runOptimizeLoop_(/*increase=*/ true); +} + + +// Called when we should use a fixed ALQ value +template +void +GasLiftSingleWell:: +updateWellStateAlqFixedValue_(const GasLiftOpt::Well &well) +{ + auto& max_alq_optional = well.max_rate(); + if (max_alq_optional) { + // According to WLIFTOPT, item 3: + // If item 2 is NO, then item 3 is regarded as the fixed + // lift gas injection rate for the well. + auto new_alq = *max_alq_optional; + this->well_state_.setALQ(this->well_name_, new_alq); + } + // else { + // // If item 3 is defaulted, the lift gas rate remains + // // unchanged at its current value. + //} + +} + +// Determine if we should use a fixed ALQ value. +// +// From the manual for WLIFTOPT, item 2: +// Is the well's lift gas injection rate to be calculated by the +// optimization facility? +// - YES : The well's lift gas injection rate is calculated by the +// optimization facility. +// - NO : The well's lift gas injection rate remains fixed at a +// value that can be set either in Item 3 of this keyword, or in +// Item 12 of keyword WCONPROD, or with keyword WELTARG. +template +bool +GasLiftSingleWell:: +useFixedAlq_(const GasLiftOpt::Well &well) +{ + auto wliftopt_item2 = well.use_glo(); + if (wliftopt_item2) { + return false; + } + else { + // auto& max_alq_optional = well.max_rate(); + // if (max_alq_optional) { + // According to WLIFTOPT, item 3: + // If item 2 is NO, then item 3 is regarded as the fixed + // lift gas injection rate for the well. + // } + // else { + // If item 3 is defaulted, the lift gas rate remains + // unchanged at its current value. + // } + return true; + } +} + +template +void +GasLiftSingleWell:: +warnMaxIterationsExceeded_() +{ + const std::string msg = fmt::format( + "Max iterations ({}) exceeded", this->max_iterations_); + displayWarning_(msg); +} + +/**************************************** + * Methods declared in OptimizeState + ****************************************/ + +template +std::pair, bool> +GasLiftSingleWell::OptimizeState:: +addOrSubtractAlqIncrement(double alq) +{ + auto [alq_opt, limited] + = this->parent.addOrSubtractAlqIncrement_(alq, this->increase); + if (!alq_opt) { + const std::string msg = fmt::format( + "iteration {}, alq = {} : not able to {} ALQ increment", + this->it, alq, (this->increase ? "add" : "subtract")); + } + return {alq_opt, limited}; +} + +template +double +GasLiftSingleWell::OptimizeState:: +calcEcoGradient( + double oil_rate, double new_oil_rate, double gas_rate, double new_gas_rate) +{ + return this->parent.calcEcoGradient_( + oil_rate, new_oil_rate, gas_rate, new_gas_rate, this->increase); +} + +// NOTE: According to WLIFTOPT item 5 : +// if min_rate() is negative, it means: allocate at least enough lift gas +// to enable the well to flow +// We will interpret this as (see discussion above GasLiftSingleWell() +// in this file): Allocate at least the amount of lift gas needed to +// get a positive oil production rate. +template +bool +GasLiftSingleWell::OptimizeState:: +checkAlqOutsideLimits(double alq, double oil_rate) +{ + std::ostringstream ss; + bool result = false; + + if (this->increase) { + if (alq >= this->parent.max_alq_) { + ss << "ALQ >= " << this->parent.max_alq_ << " (max limit), " + << "stopping iteration"; + result = true; + } + else { // checking the minimum limit... + // NOTE: A negative min_alq_ means: allocate at least enough lift gas + // to enable the well to flow, see WLIFTOPT item 5. + if (this->parent.min_alq_ < 0) { + // - if oil rate is negative (i.e. the well is not flowing), continue to + // increase ALQ (according WLIFTOPT item 5) and try make the well + // flow. + // - else if oil rate is already positive, there is no minimum + // limit for ALQ in this case + result = false; + } + else { + // NOTE: checking for a lower limit is not necessary + // when increasing alq. If ALQ was smaller than the minimum when + // we entered the runOptimizeLoop_() method, + // increaseALQtoMinALQ_() will ensure that ALQ >= min_alq + assert(alq >= this->parent.min_alq_ ); + result = false; + } + } + } + else { // we are decreasing lift gas + if ( alq == 0 ) { + ss << "ALQ is zero, cannot decrease further. Stopping iteration."; + return true; + } + else if ( alq < 0 ) { + ss << "Negative ALQ: " << alq << ". Stopping iteration."; + return true; + } + // NOTE: A negative min_alq_ means: allocate at least enough lift gas + // to enable the well to flow, see WLIFTOPT item 5. + if (this->parent.min_alq_ < 0) { + // We know that the well is flowing (oil_rate > 0) since that was + // already checked in runOptimizeLoop_() by calling checkNegativeOilRate() + assert(oil_rate >= 0); + result = false; + } + else { + if (alq <= this->parent.min_alq_ ) { + // According to WLIFTOPT item 5: + // "If a positive value is specified, the well is + // allocated at least that amount of lift gas, + // unless the well is unable to flow with that rate + // of lift gas injection, or unless the well can + // already meet one of its own rate limits before + // receiving its minimum lift gas rate." + // + // - We already know that the well is flowing, (oil_rate > 0), + // since that was already checked in runOptimizeLoop_() by calling + // checkNegativeOilRate(). + // - We also know that the rate limit was not exceeded since that was + // checked by checkWellRatesViolated() + assert( oil_rate >= 0); + ss << "ALQ <= " << this->parent.min_alq_ << " (min limit), " + << "stopping iteration"; + result = true; + } + else { + // NOTE: checking for an upper limit should not be necessary + // when decreasing alq.. so this is just to catch an + // illegal state at an early point. + if (alq >= this->parent.max_alq_) { + warn_( "unexpected: alq above upper limit when trying to " + "decrease lift gas. aborting iteration."); + result = true; + } + else { + result = false; + } + } + } + } + if (this->parent.debug) { + const std::string msg = ss.str(); + if (!msg.empty()) + this->parent.displayDebugMessage_(msg); + } + return result; +} + +template +bool +GasLiftSingleWell::OptimizeState:: +checkNegativeOilRate(double oil_rate) +{ + if (oil_rate < 0) { + const std::string msg = fmt::format( + "Negative oil rate {}. Stopping iteration", oil_rate); + this->parent.displayDebugMessage_(msg); + return true; + } + return false; +} + +// +// bool checkEcoGradient(double gradient) +// +// - Determine if the gradient has reached the limit of the economic gradient. +// +// - If we are increasing lift gas, returns true if the gradient is smaller +// than or equal to the economic gradient, +// +// - If we are decreasing lift gas, returns true if the gradient is greater +// than or equal to the economic gradient. (I.e., we assume too much lift gas +// is being used and the gradient has become too small. We try to decrease +// lift gas until the gradient increases and reaches the economic gradient..) +// +template +bool +GasLiftSingleWell::OptimizeState:: +checkEcoGradient(double gradient) +{ + std::ostringstream ss; + bool result = false; + + if (this->parent.debug) { + ss << "checking gradient: " << gradient; + } + if (this->increase) { + if (this->parent.debug) ss << " <= " << this->parent.eco_grad_ << " --> "; + if (gradient <= this->parent.eco_grad_) { + if (this->parent.debug) ss << "yes, stopping"; + result = true; + } + else { + if (this->parent.debug) ss << "no, continue"; + } + } + else { // decreasing lift gas + if (this->parent.debug) ss << " >= " << this->parent.eco_grad_ << " --> "; + if (gradient >= this->parent.eco_grad_) { + if (this->parent.debug) ss << "yes, stopping"; + result = true; + } + else { + if (this->parent.debug) ss << "no, continue"; + } + } + if (this->parent.debug) this->parent.displayDebugMessage_(ss.str()); + return result; +} + +template +bool +GasLiftSingleWell::OptimizeState:: +checkRate(double rate, double limit, const std::string &rate_str) const +{ + if (limit < rate) { + if (this->parent.debug) { + const std::string msg = fmt::format( + "iteration {} : {} rate {} exceeds target {}. Stopping iteration", + this->it, rate_str, rate, limit); + this->parent.displayDebugMessage_(msg); + } + return true; + } + return false; +} + +template +bool +GasLiftSingleWell::OptimizeState:: +checkWellRatesViolated(std::vector &potentials) +{ + auto callback = [*this](double rate, double limit, const std::string &rate_str) + -> bool + { return this->checkRate(rate, limit, rate_str); }; + return this->parent.checkWellRatesViolated_( + potentials, callback, this->increase); +} + +template +bool +GasLiftSingleWell::OptimizeState:: +computeBhpAtThpLimit(double alq) +{ + auto bhp_opt = this->parent.computeBhpAtThpLimit_(alq); + if (bhp_opt) { + this->bhp = *bhp_opt; + return true; + } + else { + return false; + } +} + +template +void +GasLiftSingleWell::OptimizeState:: +debugShowIterationInfo(double alq) +{ + const std::string msg = fmt::format("iteration {}, ALQ = {}", this->it, alq); + this->parent.displayDebugMessage_(msg); +} + + +// NOTE: When calculating the gradient, determine what the well would produce if +// the lift gas injection rate were increased by one increment. The +// production rates are adjusted if necessary to obey +// any rate or BHP limits that the well may be subject to. From this +// information, calculate the well's "weighted incremental +// gradient" +// +// TODO: What does it mean to "adjust the production rates" given a +// BHP limit? +// +template +double +GasLiftSingleWell::OptimizeState:: +getBhpWithLimit() +{ + auto [new_bhp, limited] = this->parent.getBhpWithLimit_(this->bhp); + if (limited) { + // TODO: is it possible that bhp falls below the limit when + // adding lift gas? I.e. if this->increase == true.. + // TODO: we keep the current alq, but it should probably + // be adjusted since we changed computed bhp. But how? + + // Stop iteration, but first check the economic gradient + // with the bhp_update. If the gradient looks OK (see the + // main optimize loop) we keep the current ALQ value. + this->stop_iteration = true; + } + return new_bhp; +} diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp new file mode 100644 index 000000000..aded0b167 --- /dev/null +++ b/opm/simulators/wells/GasLiftStage2.hpp @@ -0,0 +1,225 @@ +/* + Copyright 2021 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_GASLIFT_STAGE2_HEADER_INCLUDED +#define OPM_GASLIFT_STAGE2_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +// NOTE: BlackoilWellModel.hpp includes ourself (GasLiftStage2.hpp), so we need +// to forward declare BlackoilWellModel for it to be defined in this file. +namespace Opm { + template class BlackoilWellModel; +} +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + template + class GasLiftStage2 { + using Simulator = GetPropType; + using WellState = WellStateFullyImplicitBlackoil; + using BlackoilWellModel = Opm::BlackoilWellModel; + using GasLiftSingleWell = Opm::GasLiftSingleWell; + using GLiftWellState = Opm::GasLiftWellState; + using GLiftOptWells = typename BlackoilWellModel::GLiftOptWells; + using GLiftProdWells = typename BlackoilWellModel::GLiftProdWells; + using GLiftWellStateMap = typename BlackoilWellModel::GLiftWellStateMap; + using GradPair = std::pair; + using GradPairItr = std::vector::iterator; + using GradInfo = typename GasLiftSingleWell::GradInfo; + using GradMap = std::map; + static const int Water = BlackoilPhases::Aqua; + static const int Oil = BlackoilPhases::Liquid; + static const int Gas = BlackoilPhases::Vapour; + public: + GasLiftStage2( + const BlackoilWellModel &well_model, + const Simulator &ebos_simulator, + DeferredLogger &deferred_logger, + WellState &well_state, + GLiftProdWells &prod_wells, + GLiftOptWells &glift_wells, + GLiftWellStateMap &state_map + ); + void runOptimize(); + private: + void addOrRemoveALQincrement_( + GradMap &grad_map, const std::string well_name, bool add); + std::optional calcIncOrDecGrad_( + const std::string name, const GasLiftSingleWell &gs_well, bool increase); + bool checkRateAlreadyLimited_(GLiftWellState &state, bool increase); + GradInfo deleteDecGradItem_(const std::string &name); + GradInfo deleteIncGradItem_(const std::string &name); + GradInfo deleteGrad_(const std::string &name, bool increase); + void displayDebugMessage_(const std::string &msg); + void displayDebugMessage2B_(const std::string &msg); + void displayDebugMessage_(const std::string &msg, const std::string &group_name); + void displayWarning_(const std::string &msg, const std::string &group_name); + void displayWarning_(const std::string &msg); + std::tuple getCurrentGroupRates_( + const Opm::Group &group); + std::tuple getCurrentGroupRatesRecursive_( + const Opm::Group &group); + std::tuple getCurrentWellRates_( + const std::string &well_name, const std::string &group_name); + std::vector getGroupGliftWells_( + const Opm::Group &group); + void getGroupGliftWellsRecursive_( + const Opm::Group &group, std::vector &wells); + std::pair getStdWellRates_(const WellInterface &well); + void optimizeGroup_(const Opm::Group &group); + void optimizeGroupsRecursive_(const Opm::Group &group); + void recalculateGradientAndUpdateData_( + GradPairItr &grad_itr, bool increase, + std::vector &grads, std::vector &other_grads); + void redistributeALQ_( + std::vector &wells, const Opm::Group &group, + std::vector &inc_grads, std::vector &dec_grads); + void removeSurplusALQ_( + const Opm::Group &group, + std::vector &inc_grads, std::vector &dec_grads); + void saveGrad_(GradMap &map, const std::string &name, GradInfo &grad); + void saveDecGrad_(const std::string &name, GradInfo &grad); + void saveIncGrad_(const std::string &name, GradInfo &grad); + void sortGradients_(std::vector &grads); + std::optional updateGrad_( + const std::string &name, GradInfo &grad, bool increase); + void updateGradVector_( + const std::string &name, std::vector &grads, double grad); + + DeferredLogger &deferred_logger_; + const Simulator &ebos_simulator_; + const BlackoilWellModel &well_model_; + WellState &well_state_; + GLiftProdWells &prod_wells_; + GLiftOptWells &stage1_wells_; + GLiftWellStateMap &well_state_map_; + + int report_step_idx_; + const SummaryState &summary_state_; + const Schedule &schedule_; + const PhaseUsage &phase_usage_; + const GasLiftOpt& glo_; + GradMap inc_grads_; + GradMap dec_grads_; + bool debug_; + int max_iterations_ = 1000; + //int time_step_idx_; + int nonlinear_iteration_idx_; + + struct OptimizeState { + OptimizeState( GasLiftStage2 &parent_, const Opm::Group &group_ ) : + parent{parent_}, + group{group_}, + it{0} + {} + GasLiftStage2 &parent; + const Opm::Group &group; + int it; + + using GradInfo = typename GasLiftStage2::GradInfo; + using GradPair = typename GasLiftStage2::GradPair; + using GradPairItr = typename GasLiftStage2::GradPairItr; + using GradMap = typename GasLiftStage2::GradMap; + void calculateEcoGradients(std::vector &wells, + std::vector &inc_grads, std::vector &dec_grads); + bool checkAtLeastTwoWells(std::vector &wells); + void debugShowIterationInfo(); + std::pair,std::optional> + getEcoGradients( + std::vector &inc_grads, std::vector &dec_grads); + void recalculateGradients( + std::vector &inc_grads, std::vector &dec_grads, + GradPairItr &min_dec_grad_itr, GradPairItr &max_inc_grad_itr); + void redistributeALQ( GradPairItr &min_dec_grad, GradPairItr &max_inc_grad); + + private: + void displayDebugMessage_(const std::string &msg); + void displayWarning_(const std::string &msg); + + }; + + struct SurplusState { + SurplusState( GasLiftStage2 &parent_, const Opm::Group &group_, + double oil_rate_, double gas_rate_, double alq_, double min_eco_grad_, + double oil_target_, double gas_target_, + std::optional max_glift_) : + parent{parent_}, + group{group_}, + oil_rate{oil_rate_}, + gas_rate{gas_rate_}, + alq{alq_}, + min_eco_grad{min_eco_grad_}, + oil_target{oil_target_}, + gas_target{gas_target_}, + max_glift{max_glift_}, + it{0} + {} + GasLiftStage2 &parent; + const Opm::Group &group; + double oil_rate; + double gas_rate; + double alq; + const double min_eco_grad; + const double oil_target; + const double gas_target; + std::optional max_glift; + int it; + + void addOrRemoveALQincrement( + GradMap &grad_map, const std::string well_name, bool add); + bool checkALQlimit(); + bool checkEcoGradient(const std::string &well_name, double eco_grad); + bool checkGasTarget(); + bool checkOilTarget(); + void updateRates(const std::string &name, const GradInfo &gi); + private: + }; + }; + +#include "GasLiftStage2_impl.hpp" + +} // namespace Opm + +#endif // OPM_GASLIFT_STAGE2_HEADER_INCLUDED diff --git a/opm/simulators/wells/GasLiftStage2_impl.hpp b/opm/simulators/wells/GasLiftStage2_impl.hpp new file mode 100644 index 000000000..14026882d --- /dev/null +++ b/opm/simulators/wells/GasLiftStage2_impl.hpp @@ -0,0 +1,1027 @@ +/* + Copyright 2021 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +template +GasLiftStage2:: +GasLiftStage2( + const BlackoilWellModel &well_model, + const Simulator &ebos_simulator, + DeferredLogger &deferred_logger, + WellState &well_state, + GLiftProdWells &prod_wells, + GLiftOptWells &glift_wells, + GLiftWellStateMap &state_map +) : + deferred_logger_{deferred_logger}, + ebos_simulator_{ebos_simulator}, + well_model_{well_model}, + well_state_{well_state}, + prod_wells_{prod_wells}, + stage1_wells_{glift_wells}, + well_state_map_{state_map}, + report_step_idx_{ebos_simulator_.episodeIndex()}, + summary_state_{ebos_simulator_.vanguard().summaryState()}, + schedule_{ebos_simulator.vanguard().schedule()}, + phase_usage_{well_model_.phaseUsage()}, + glo_{schedule_.glo(report_step_idx_)}, + debug_{false} +{ +// this->time_step_idx_ +// = this->ebos_simulator_.model().newtonMethod().currentTimeStep(); + this->nonlinear_iteration_idx_ + = this->ebos_simulator_.model().newtonMethod().numIterations(); + +} + +/******************************************** + * Public methods in alphabetical order + ********************************************/ + +// runOptimize(): +// +// If a group has any production rate constraints, and/or a limit on +// its total rate of lift gas supply, allocates lift gas +// preferentially to the wells that gain the most benefit from +// it. Lift gas increments are allocated in turn to the well that +// currently has the largest weighted incremental gradient. The +// procedure takes account of any limits on the group production rate +// or lift gas supply applied to any level of group, including the FIELD level group. +template +void +GasLiftStage2:: +runOptimize() +{ + const auto& group = this->schedule_.getGroup("FIELD", this->report_step_idx_); + + optimizeGroupsRecursive_(group); + +} + + +/******************************************** + * Private methods in alphabetical order + ********************************************/ + +// 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") +template +void +GasLiftStage2:: +addOrRemoveALQincrement_(GradMap &grad_map, const std::string well_name, bool add) +{ + const GradInfo &gi = grad_map.at(well_name); + GLiftWellState &state = *(this->well_state_map_.at(well_name).get()); + if (this->debug_) { + auto new_alq = gi.alq; + auto old_alq = state.alq(); + const std::string msg = fmt::format("well {} : {} ALQ increment, " + "old alq: {}, new alq: {}", + 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, add); + //gs_well.updateStage2State(gi, add); + this->well_state_.setALQ(well_name, gi.alq); +} + +template +std::optional::GradInfo> +GasLiftStage2:: +calcIncOrDecGrad_( + const std::string well_name, const GasLiftSingleWell &gs_well, bool increase) +{ + GLiftWellState &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); + */ + return std::nullopt; + } + else { + auto [oil_rate, gas_rate] = state.getRates(); + auto alq = state.alq(); + auto grad = gs_well.calcIncOrDecGradient(oil_rate, gas_rate, alq, increase); + if (grad) { + const std::string msg = fmt::format( + "well {} : adding {} gradient = {}", + well_name, + (increase ? "incremental" : "decremental"), + grad->grad + ); + displayDebugMessage_(msg); + } + return grad; + } +} + +template +bool +GasLiftStage2:: +checkRateAlreadyLimited_(GLiftWellState &state, bool increase) +{ + auto current_increase = state.increase(); + bool do_check = false; + if (current_increase) { + if (*current_increase == increase) do_check = true; + } + else { + // If current_increase is not defined, it means that stage1 + // was unable to either increase nor decrease the ALQ. If the + // initial rates stored in "state" is limited, and if + // "increase" is true, it is not likely that adding ALQ will + // cause the new rates not to be limited. However, if + // "increase" is false, subtracting ALQ can make the new rates + // not limited. + if (increase) do_check = true; + } + if (do_check) { + if (state.gasIsLimited() || state.oilIsLimited() || state.alqIsLimited()) { + const std::string msg = fmt::format( + "{} gradient : skipping since {} was limited in previous step", + (increase ? "incremental" : "decremental"), + (state.oilIsLimited() ? "oil" : + (state.gasIsLimited() ? "gas" : "alq"))); + displayDebugMessage_(msg); + return true; + } + } + return false; +} + + +template +typename GasLiftStage2::GradInfo +GasLiftStage2:: +deleteGrad_(const std::string &name, bool increase) +{ + GradMap &map = increase ? this->inc_grads_ : this->dec_grads_; + auto value = map.at(name); + map.erase(name); + return value; +} + +template +typename GasLiftStage2::GradInfo +GasLiftStage2:: +deleteDecGradItem_(const std::string &name) +{ + return deleteGrad_(name, /*increase=*/false); +} + +template +typename GasLiftStage2::GradInfo +GasLiftStage2:: +deleteIncGradItem_(const std::string &name) +{ + return deleteGrad_(name, /*increase=*/true); +} + +template +void +GasLiftStage2:: +displayWarning_(const std::string &msg, const std::string &group_name) +{ + const std::string message = fmt::format( + "GAS LIFT OPTIMIZATION (STAGE2), GROUP: {} : {}", group_name, msg); + this->deferred_logger_.warning("WARNING", message); +} + +template +void +GasLiftStage2:: +displayWarning_(const std::string &msg) +{ + const std::string message = fmt::format( + "GAS LIFT OPTIMIZATION (STAGE2) : {}", msg); + this->deferred_logger_.warning("WARNING", message); +} + +template +void +GasLiftStage2:: +displayDebugMessage_(const std::string &msg) +{ + if (this->debug_) { + const std::string message = fmt::format( + " GLIFT2 (DEBUG) : {}", msg); + this->deferred_logger_.info(message); + } +} + +template +void +GasLiftStage2:: +displayDebugMessage2B_(const std::string &msg) +{ + if (this->debug_) { + const std::string message = fmt::format( + "Stage 2B : {}", msg); + displayDebugMessage_(message); + } +} +template +void +GasLiftStage2:: +displayDebugMessage_(const std::string &msg, const std::string &group_name) +{ + if (this->debug_) { + const std::string message = fmt::format( + "Group {} : {}", group_name, msg); + displayDebugMessage_(message); + } +} + +template +std::tuple +GasLiftStage2:: +getCurrentGroupRates_(const Opm::Group &group) +{ + if (this->debug_) { + auto [oil_rate, gas_rate, alq] = getCurrentGroupRatesRecursive_(group); + const std::string msg = fmt::format( + "Current group rates for {} : oil: {}, gas: {}, alq: {}", + group.name(), oil_rate, gas_rate, alq); + displayDebugMessage2B_(msg); + return {oil_rate, gas_rate, alq}; + } + else { + return getCurrentGroupRatesRecursive_(group); + } +} + + +template +std::tuple +GasLiftStage2:: +getCurrentGroupRatesRecursive_(const Opm::Group &group) +{ + double oil_rate = 0.0; + double gas_rate = 0.0; + double alq = 0.0; + // NOTE: A group can either contain wells or groups, but not both + if (group.wellgroup()) { + for (const std::string& well_name : group.wells()) { + auto [sw_oil_rate, sw_gas_rate, sw_alq] = + getCurrentWellRates_(well_name, group.name()); + oil_rate += sw_oil_rate; + gas_rate += sw_gas_rate; + alq += sw_alq; + } + } + else { + for (const std::string& group_name : group.groups()) { + if(this->schedule_.back().groups.has(group_name)) { + const Group& sub_group = + this->schedule_.getGroup(group_name, this->report_step_idx_); + // If groups have efficiency factors to model + // synchronized downtime of their subordinate wells + // (see keyword GEFAC), their lift gas injection rates + // are multiplied by their efficiency factors when + // they are added to the lift gas supply rate of the + // parent group. + const auto gefac = sub_group.getGroupEfficiencyFactor(); + auto [sg_oil_rate, sg_gas_rate, sg_alq] = + getCurrentGroupRatesRecursive_(sub_group); + oil_rate += (gefac * sg_oil_rate); + gas_rate += (gefac * sg_gas_rate); + alq += (gefac * sg_alq); + } + } + } + return {oil_rate, gas_rate, alq}; +} + +template +std::tuple +GasLiftStage2:: +getCurrentWellRates_(const std::string &well_name, const std::string &group_name) +{ + double oil_rate, gas_rate, alq; + bool success = false; + const WellInterface *well_ptr = nullptr; + std::string debug_info; + if (this->stage1_wells_.count(well_name) == 1) { + GasLiftSingleWell &gs_well = *(this->stage1_wells_.at(well_name).get()); + const WellInterface &well = gs_well.getStdWell(); + well_ptr = &well; + GLiftWellState &state = *(this->well_state_map_.at(well_name).get()); + std::tie(oil_rate, gas_rate) = state.getRates(); + success = true; + if ( this->debug_) debug_info = "(A)"; + } + else if (this->prod_wells_.count(well_name) == 1) { + well_ptr = this->prod_wells_.at(well_name); + std::tie(oil_rate, gas_rate) = getStdWellRates_(*well_ptr); + success = true; + if ( this->debug_) debug_info = "(B)"; + } + if (success) { + assert(well_ptr); + assert(well_ptr->isProducer()); + alq = this->well_state_.getALQ(well_name); + if (this->debug_) { + const std::string msg = fmt::format( + "Rates {} for well {} : oil: {}, gas: {}, alq: {}", + debug_info, well_name, oil_rate, gas_rate, alq); + displayDebugMessage_(msg, group_name); + } + // If wells have efficiency factors to take account of regular + // downtime (see keyword WEFAC), their lift gas injection + // rates are multiplied by their efficiency factors when they + // are added to the group lift gas supply rate. This is + // consistent with the summation of flow rates for wells with + // downtime, and preserves the ratio of production rate to + // lift gas injection rate. + const auto &well_ecl = well_ptr->wellEcl(); + double factor = well_ecl.getEfficiencyFactor(); + oil_rate *= factor; + gas_rate *= factor; + alq *= factor; + if (this->debug_ && (factor != 1)) { + const std::string msg = fmt::format( + "Well {} : efficiency factor {}. New rates : oil: {}, gas: {}, alq: {}", + well_name, factor, oil_rate, gas_rate, alq); + displayDebugMessage_(msg, group_name); + } + } + else { + // NOTE: This happens for wells that are not producers, or not active. + if (this->debug_) { + const std::string msg = fmt::format("Could not determine current rates for " + "well {}: (not active or injector)", well_name); + displayDebugMessage_(msg, group_name); + } + oil_rate = 0.0; gas_rate = 0.0; alq = 0.0; + } + return std::make_tuple(oil_rate, gas_rate, alq); +} + +template +std::pair +GasLiftStage2:: +getStdWellRates_(const WellInterface &well) +{ + const int well_index = well.indexOfWell(); + const int np = this->well_state_.numPhases(); + const auto& pu = well.phaseUsage(); + auto oil_rate = + -this->well_state_.wellRates()[np * well_index + pu.phase_pos[Oil]]; + auto gas_rate = + -this->well_state_.wellRates()[np * well_index + pu.phase_pos[Gas]]; + return {oil_rate, gas_rate}; +} + +// Find all subordinate wells of a given group. +// +// NOTE: A group can either contain wells or groups, not both. +// If it contains groups, we have to traverse those recursively to find the wells. +// +// NOTE: This means that wells are located at the leaf nodes of the tree, and +// groups are located at the other nodes (not leaf nodes) of the tree +// +template +std::vector *> +GasLiftStage2:: +getGroupGliftWells_(const Opm::Group &group) +{ + std::vector wells; + getGroupGliftWellsRecursive_(group, wells); + return wells; +} + +template +void +GasLiftStage2:: +getGroupGliftWellsRecursive_(const Opm::Group &group, + std::vector &wells) +{ + for (const std::string& group_name : group.groups()) { + if(this->schedule_.back().groups.has(group_name)) { + const Group& sub_group = + this->schedule_.getGroup(group_name, this->report_step_idx_); + getGroupGliftWellsRecursive_(sub_group, wells); + } + } + for (const std::string& well_name : group.wells()) { + if (this->stage1_wells_.count(well_name) == 1) { + GasLiftSingleWell *well_ptr = this->stage1_wells_.at(well_name).get(); + wells.push_back(well_ptr); + } + } +} + +template +void +GasLiftStage2:: +optimizeGroup_(const Opm::Group &group) +{ + const auto &gl_group = this->glo_.group(group.name()); + const auto &max_glift = gl_group.max_lift_gas(); + const auto &max_total_gas = gl_group.max_total_gas(); + if (group.has_control(Group::ProductionCMode::ORAT) + || max_glift || max_total_gas) + { + displayDebugMessage_("optimizing", group.name()); + auto wells = getGroupGliftWells_(group); + std::vector inc_grads; + std::vector dec_grads; + redistributeALQ_(wells, group, inc_grads, dec_grads); + removeSurplusALQ_(group, inc_grads, dec_grads); + } + else { + displayDebugMessage_("skipping", group.name()); + } +} + +template +void +GasLiftStage2:: +optimizeGroupsRecursive_(const Opm::Group &group) +{ + for (const std::string& group_name : group.groups()) { + if(!this->schedule_.back().groups.has(group_name)) + continue; + const Group& sub_group = this->schedule_.getGroup( + group_name, this->report_step_idx_); + optimizeGroupsRecursive_(sub_group); + } + // TODO: should we also optimize groups that do not have GLIFTOPT defined? + // (i.e. glo_.has_group(name) returns false) + // IF GLIFTOPT is not specified for the group or if item 2 of GLIFTOPT + // is defaulted, there is no maximum lift gas supply for the group. + // But even if there is no limit on the liftgas supply it can still + // be desireable to use as little ALQ as possible to achieve a + // group oil rate limit or gas rate limit. + if (this->glo_.has_group(group.name())) // only optimize if GLIFTOPT is given + optimizeGroup_(group); + +} + +template +void +GasLiftStage2:: +recalculateGradientAndUpdateData_( + GradPairItr &grad_itr, bool increase, + 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) + const std::string name = grad_itr->first; + GasLiftSingleWell &gs_well = *(this->stage1_wells_.at(name).get()); + auto grad = calcIncOrDecGrad_(name, gs_well, increase); + std::optional old_grad = std::nullopt; + 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); + } + 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); + updateGradVector_(name, other_grads, old_grad->grad); + } +} + +// redistributeALQ_() : +// +// DESCRIPTION: The currently allocated increments are redistributed by +// comparing the largest weighted incremental gradient and the +// smallest weighted decremental gradient of the subordinate +// wells. Consider, for example, that the largest weighted incremental +// gradient belongs to well W1, and exceeds the smallest weighted +// decremental gradient, which belongs to well W2. It would be more +// profitable to subtract an increment of lift gas from well W2 and +// allocate it to well W1. After the exchange, the incremental and +// decremental gradients of wells W1 and W2 are recalculated. The +// exchange of increments continues until no weighted incremental +// gradient exceeds a weighted decremental gradient. +// +// TODO: maybe the redistribution can be simplified by only doing it for the +// FIELD group: For example: +// +// +// +// FIELD +// | +// PLAT-A +// ---------------+--------------------- +// | | +// M5S M5N +// ---------+---------- -----+------- +// | | | | +// B1 G1 C1 F1 +// ----+------ ---+--- ---+--- ---+--- +// | | | | | | | | | +// B-1H B-2H B-3H G-3H G-4H C-1H C-2H F-1H F-2H +// +// +// it is probably unecessary to first redistribute ALQ for the wells B-1H, B-2H, +// and B-3H in group B1, and then in a next separate step, redistribute ALQ for the +// wells G-3H, and G-4H, and similarly, for the wells in groups C1, and F1, +// and then, for the wells under M5S, and then M5N, and finally repeat the procedure +// for all the wells under group PLAT-A, i.e. the wells: +// B-1H, B-2H, B-3H, G-3H, G-4H, C-1H, C-2H, F-1H, and F-2H. +// It seems like it would be more efficient to +// just do it once for the topmost group "PLAT-A" and then skip redistribution for +// all sub groups of "PLAT-A" +// +template +void +GasLiftStage2:: +redistributeALQ_(std::vector &wells, const Opm::Group &group, + std::vector &inc_grads, std::vector &dec_grads) +{ + OptimizeState state {*this, group}; + inc_grads.reserve(wells.size()); + dec_grads.reserve(wells.size()); + state.calculateEcoGradients(wells, inc_grads, dec_grads); + if (!state.checkAtLeastTwoWells(wells)) { + // NOTE: Even though we here in redistributeALQ_() do not use the + // economic gradient if there is only a single well, we still + // need to calculate it since inc_grads and dec_grads are returned + // and will be used by removeSurplusALQ_() later. + return; + } + assert(wells.size() >= 2); + bool stop_iteration = false; + while (!stop_iteration && (state.it++ <= this->max_iterations_)) { + state.debugShowIterationInfo(); + auto [min_dec_grad, max_inc_grad] + = state.getEcoGradients(inc_grads, dec_grads); + if (min_dec_grad) { + assert( max_inc_grad ); + // Redistribute if the largest incremental gradient exceeds the + // smallest decremental gradient + if ((*max_inc_grad)->second > (*min_dec_grad)->second) { + state.redistributeALQ(*min_dec_grad, *max_inc_grad); + state.recalculateGradients( + inc_grads, dec_grads, *min_dec_grad, *max_inc_grad); + continue; + } + } + stop_iteration = true; + } + if (state.it > this->max_iterations_) { + const std::string msg = fmt::format("Max iterations {} exceeded.", + this->max_iterations_); + displayWarning_(msg, group.name()); + } +} + +// The group has surplus lift gas if it exceeds any production rate limits +// or a lift gas supply limit, or contains any wells that have a weighted +// decremental gradient less than the minimum economic gradient. +// Lift gas increments are removed in turn from the well that currently has +// the smallest weighted decremental gradient, until there is no surplus +// lift gas in the group. +template +void +GasLiftStage2:: +removeSurplusALQ_(const Opm::Group &group, + std::vector &inc_grads, std::vector &dec_grads) +{ + if (dec_grads.size() == 0) { + displayDebugMessage2B_("no wells to remove ALQ from. Skipping"); + return; + } + assert(dec_grads.size() > 0); + const auto &gl_group = this->glo_.group(group.name()); + const auto &max_glift = gl_group.max_lift_gas(); + const auto controls = group.productionControls(this->summary_state_); + //const auto &max_total_gas = gl_group.max_total_gas(); + auto [oil_rate, gas_rate, alq] = getCurrentGroupRates_(group); + auto min_eco_grad = this->glo_.min_eco_gradient(); + bool stop_iteration = false; + 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: {}. " + "oil_rate = {}, oil_target = {}, gas_rate = {}, gas_target = {}, " + "alq = {}, max_alq = {}", group.name(), oil_rate, controls.oil_target, + gas_rate, controls.gas_target, alq, max_glift_str); + displayDebugMessage2B_(msg); + } + SurplusState state {*this, group, oil_rate, gas_rate, alq, + min_eco_grad, controls.oil_target, controls.gas_target, max_glift }; + while (!stop_iteration) { + if (dec_grads.size() >= 2) { + sortGradients_(dec_grads); + } + auto dec_grad_itr = dec_grads.begin(); + const auto &well_name = dec_grad_itr->first; + auto eco_grad = dec_grad_itr->second; + bool remove = false; + if (state.checkOilTarget() || state.checkGasTarget() || state.checkALQlimit()) { + remove = true; + } + else { + // NOTE: It is enough to check the economic gradient of the first well + // in dec_grads since they are sorted according to the eco. grad. + // If the first well's eco. grad. is greater than the minimum + // eco. grad. then all the other wells' eco. grad. will also be + // greater. + if (state.checkEcoGradient(well_name, eco_grad)) remove = true; + } + if (remove) { + const GradInfo &gi = this->dec_grads_.at(well_name); + state.updateRates(well_name, gi); + state.addOrRemoveALQincrement( this->dec_grads_, well_name, /*add=*/false); + recalculateGradientAndUpdateData_( + dec_grad_itr, /*increase=*/false, dec_grads, 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 + if (dec_grads.size() == 0) stop_iteration = true; + ++state.it; + } + else { + stop_iteration = true; + } + } + if (state.it >= 1) { + if (this->debug_) { + auto [oil_rate2, gas_rate2, alq2] = getCurrentGroupRates_(group); + const std::string msg = fmt::format( + "Finished after {} iterations for group: {}." + " oil_rate = {}, gas_rate = {}, alq = {}", state.it, + group.name(), oil_rate2, gas_rate2, alq2); + displayDebugMessage2B_(msg); + } + } + else { + displayDebugMessage2B_("Finished after 0 iterations"); + } +} + +template +void +GasLiftStage2:: +saveGrad_(GradMap &map, const std::string &name, GradInfo &grad) +{ + if (auto it = map.find(name); it == map.end()) { + auto [map_it, success] = map.emplace(name, grad); + assert(success); + } + else { + it->second = grad; + } +} + +template +void +GasLiftStage2:: +saveDecGrad_(const std::string &name, GradInfo &grad) +{ + saveGrad_(this->dec_grads_, name, grad); +} + +template +void +GasLiftStage2:: +saveIncGrad_(const std::string &name, GradInfo &grad) +{ + saveGrad_(this->inc_grads_, name, grad); +} + +template +void +GasLiftStage2:: +sortGradients_(std::vector &grads) +{ + auto cmp = [](GradPair a, GradPair b) { + return a.second < b.second; + }; + std::sort(grads.begin(), grads.end(), cmp); +} + +template +std::optional::GradInfo> +GasLiftStage2:: +updateGrad_(const std::string &name, GradInfo &grad, bool increase) +{ + GradMap &map = increase ? this->inc_grads_ : this->dec_grads_; + std::optional old_value = std::nullopt; + if (map.count(name) == 1) { + old_value = map.at(name); + } + map[name] = grad; + return old_value; +} + +template +void +GasLiftStage2:: +updateGradVector_(const std::string &name, std::vector &grads, double grad) +{ + for (auto itr = grads.begin(); itr != grads.end(); itr++) { + if (itr->first == name) { + itr->second = grad; + return; + } + } + grads.push_back({name, grad}); + // NOTE: the gradient vector is no longer sorted, but sorting will be done + // later in getEcoGradients() +} + + +/*********************************************** + * Public methods declared in OptimizeState + ***********************************************/ + +template +void +GasLiftStage2::OptimizeState:: +calculateEcoGradients(std::vector &wells, + std::vector &inc_grads, std::vector &dec_grads) +{ + for (auto well_ptr : wells) { + const auto &gs_well = *well_ptr; // gs = GasLiftSingleWell + const auto &name = gs_well.name(); + auto inc_grad = this->parent.calcIncOrDecGrad_(name, gs_well, /*increase=*/true); + if (inc_grad) { + inc_grads.emplace_back(std::make_pair(name, inc_grad->grad)); + this->parent.saveIncGrad_(name, *inc_grad); + } + auto dec_grad = this->parent.calcIncOrDecGrad_(name, gs_well, /*increase=*/false); + if (dec_grad) { + dec_grads.emplace_back(std::make_pair(name, dec_grad->grad)); + this->parent.saveDecGrad_(name, *dec_grad); + } + } +} + + +template +bool +GasLiftStage2::OptimizeState:: +checkAtLeastTwoWells(std::vector &wells) +{ + if (wells.size() < 2) { + const std::string msg = fmt::format( + "skipping: too few wells ({}) to optimize.", wells.size()); + displayDebugMessage_(msg); + return false; + } + return true; +} + +template +void +GasLiftStage2::OptimizeState:: +debugShowIterationInfo() +{ + const std::string msg = fmt::format("iteration {}", this->it); + displayDebugMessage_(msg); +} + +template +std::pair::GradPairItr>, + std::optional::GradPairItr>> +GasLiftStage2::OptimizeState:: +getEcoGradients(std::vector &inc_grads, std::vector &dec_grads) +{ + if (inc_grads.size() > 0 && dec_grads.size() > 0) { + this->parent.sortGradients_(inc_grads); + this->parent.sortGradients_(dec_grads); + // The largest incremental gradient is the last element + auto inc_grad = std::prev(inc_grads.end()); + std::optional inc_grad_opt; + std::optional dec_grad_opt; + // The smallest decremental gradient is at the beginning + for (auto itr = dec_grads.begin(); itr != dec_grads.end(); itr++) { + if (itr->first == inc_grad->first) { + // Don't consider decremental gradients with the same well name + continue; + } + dec_grad_opt = itr; + break; + } + if (dec_grad_opt) { + inc_grad_opt = inc_grad; + return { dec_grad_opt, inc_grad_opt }; + } + } + return {std::nullopt, std::nullopt}; +} + +// Recalculate gradients (and related information, see struct GradInfo in +// GasLiftSingleWell.hpp) after an ALQ increment +// has been given from the well with minumum decremental gradient (represented +// by the input argument min_dec_grad_itr) to the well with the largest +// incremental gradient (represented by input argument max_inc_grad_itr). +// +// For the well with the largest incremental gradient, we compute a new +// incremental gradient given the new ALQ. The new decremental gradient for this +// well is set equal to the current incremental gradient (before the ALQ is added) +// Similiarly, for the well with the smalles decremental gradient, we compute +// a new decremental gradient given the new ALQ. The new incremental gradient +// for this well is set equal to the current decremental gradient +// (before the ALQ is subtracted) +template +void +GasLiftStage2::OptimizeState:: +recalculateGradients( + std::vector &inc_grads, std::vector &dec_grads, + GradPairItr &min_dec_grad_itr, GradPairItr &max_inc_grad_itr) +{ + this->parent.recalculateGradientAndUpdateData_( + max_inc_grad_itr, /*increase=*/true, inc_grads, dec_grads); + this->parent.recalculateGradientAndUpdateData_( + min_dec_grad_itr, /*increase=*/false, dec_grads, inc_grads); +} + +// Take one ALQ increment from well1, and give it to well2 +template +void +GasLiftStage2::OptimizeState:: + redistributeALQ( GradPairItr &min_dec_grad, GradPairItr &max_inc_grad) +{ + const std::string msg = fmt::format( + "redistributing ALQ from well {} (dec gradient: {}) " + "to well {} (inc gradient {})", + min_dec_grad->first, min_dec_grad->second, + max_inc_grad->first, max_inc_grad->second); + displayDebugMessage_(msg); + this->parent.addOrRemoveALQincrement_( + 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); +} + +/********************************************** + * Private methods declared in OptimizeState + **********************************************/ + +template +void +GasLiftStage2::OptimizeState:: +displayDebugMessage_(const std::string &msg) +{ + this->parent.displayDebugMessage_(msg, this->group.name()); +} + +template +void +GasLiftStage2::OptimizeState:: +displayWarning_(const std::string &msg) +{ + this->parent.displayWarning_(msg, this->group.name()); +} + +/********************************************** + * Public methods declared in SurplusState + **********************************************/ + +template +void +GasLiftStage2::SurplusState:: +addOrRemoveALQincrement(GradMap &grad_map, const std::string well_name, bool add) +{ + if (this->parent.debug_) { + const std::string msg = fmt::format("group: {} : well {} : {} ALQ increment", + this->group.name(), well_name, (add ? "adding" : "subtracting")); + this->parent.displayDebugMessage2B_(msg); + } + this->parent.addOrRemoveALQincrement_(grad_map, well_name, add); +} + +template +bool +GasLiftStage2::SurplusState:: +checkALQlimit() +{ + if (this->max_glift) { + double max_alq = *(this->max_glift); + double increment = this->parent.glo_.gaslift_increment(); + double epsilon = 1e-6 * increment; + if ((max_alq+epsilon) < this->alq ) { + if (this->parent.debug_) { + const std::string msg = fmt::format("group: {} : " + "ALQ rate {} is greater than ALQ limit {}", this->group.name(), + this->alq, max_alq); + this->parent.displayDebugMessage2B_(msg); + } + return true; + } + } + return false; +} + +template +bool +GasLiftStage2::SurplusState:: +checkEcoGradient(const std::string &well_name, double eco_grad) +{ + if (eco_grad < this->min_eco_grad) { + if (this->parent.debug_) { + const std::string msg = fmt::format("group: {}, well: {} : " + "economic gradient {} less than minimum ({})", this->group.name(), + well_name, eco_grad, this->min_eco_grad); + this->parent.displayDebugMessage2B_(msg); + } + return true; + } + else { + return false; + } +} + +template +bool +GasLiftStage2::SurplusState:: +checkGasTarget() +{ + if (this->group.has_control(Group::ProductionCMode::GRAT)) { + if (this->gas_target < this->gas_rate ) { + if (this->parent.debug_) { + const std::string msg = fmt::format("group: {} : " + "gas rate {} is greater than gas target {}", this->group.name(), + this->gas_rate, this->gas_target); + this->parent.displayDebugMessage2B_(msg); + } + return true; + } + } + return false; +} + +template +bool +GasLiftStage2::SurplusState:: +checkOilTarget() +{ + if (this->group.has_control(Group::ProductionCMode::ORAT)) { + if (this->oil_target < this->oil_rate ) { + if (this->parent.debug_) { + const std::string msg = fmt::format("group: {} : " + "oil rate {} is greater than oil target {}", this->group.name(), + this->oil_rate, this->oil_target); + this->parent.displayDebugMessage2B_(msg); + } + return true; + } + } + return false; +} + +template +void +GasLiftStage2::SurplusState:: +updateRates(const std::string &well_name, const GradInfo &gi) +{ + GLiftWellState &state = *(this->parent.well_state_map_.at(well_name).get()); + GasLiftSingleWell &gs_well = *(this->parent.stage1_wells_.at(well_name).get()); + const WellInterface &well = gs_well.getStdWell(); + const auto &well_ecl = well.wellEcl(); + double factor = well_ecl.getEfficiencyFactor(); + auto delta_oil = gi.new_oil_rate - state.oilRate(); + this->oil_rate += factor * delta_oil; + auto delta_gas = gi.new_gas_rate - state.gasRate(); + this->gas_rate += factor * delta_gas; + auto delta_alq = gi.alq - state.alq(); + this->alq += factor * delta_alq; +} diff --git a/opm/simulators/wells/GasLiftWellState.hpp b/opm/simulators/wells/GasLiftWellState.hpp new file mode 100644 index 000000000..1187b5b94 --- /dev/null +++ b/opm/simulators/wells/GasLiftWellState.hpp @@ -0,0 +1,82 @@ +/* + Copyright 2021 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + +#ifndef OPM_GASLIFT_WELL_STATE_HEADER_INCLUDED +#define OPM_GASLIFT_WELL_STATE_HEADER_INCLUDED + +#include +#include +#include +#include + +namespace Opm +{ + template + class GasLiftWellState + { + public: + GasLiftWellState() { } + GasLiftWellState(double oil_rate, bool oil_is_limited, + double gas_rate, bool gas_is_limited, + double alq, bool alq_is_limited, std::optional increase) : + oil_rate_{oil_rate}, + oil_is_limited_{oil_is_limited}, + gas_rate_{gas_rate}, + gas_is_limited_{gas_is_limited}, + alq_{alq}, + alq_is_limited_{alq_is_limited}, + increase_{increase} + {} + double alq() const { return alq_; } + bool alqChanged() { return increase_.has_value(); } + bool alqIsLimited() const { return alq_is_limited_; } + bool gasIsLimited() const { return gas_is_limited_; } + double gasRate() const { return gas_rate_; } + std::pair getRates() { return {oil_rate_, gas_rate_}; } + std::optional increase() const { return increase_; } + bool oilIsLimited() const { return oil_is_limited_; } + double oilRate() const { return oil_rate_; } + void update(double oil_rate, bool oil_is_limited, + double gas_rate, bool gas_is_limited, + double alq, bool alq_is_limited, + bool increase) + { + oil_rate_ = oil_rate; + oil_is_limited_ = oil_is_limited; + gas_rate_ = gas_rate; + gas_is_limited_ = gas_is_limited; + alq_ = alq; + alq_is_limited_ = alq_is_limited; + increase_ = increase; + } + private: + double oil_rate_; + bool oil_is_limited_; + double gas_rate_; + bool gas_is_limited_; + double alq_; + bool alq_is_limited_; + std::optional increase_; + }; + +#include "GasLiftWellState_impl.hpp" + +} // namespace Opm + +#endif // OPM_GASLIFT_WELL_STATE_HEADER_INCLUDED diff --git a/opm/simulators/wells/GasLiftWellState_impl.hpp b/opm/simulators/wells/GasLiftWellState_impl.hpp new file mode 100644 index 000000000..1a046e6d4 --- /dev/null +++ b/opm/simulators/wells/GasLiftWellState_impl.hpp @@ -0,0 +1,18 @@ +/* + Copyright 2021 Equinor ASA. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 8e192ffe5..f7e9da9c2 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -45,6 +45,10 @@ namespace Opm using typename Base::RateConverterType; using typename Base::SparseMatrixAdapter; using typename Base::FluidState; + using typename Base::GasLiftSingleWell; + using typename Base::GLiftProdWells; + using typename Base::GLiftOptWells; + using typename Base::GLiftWellStateMap; /// the number of reservior equations using Base::numEq; @@ -120,10 +124,13 @@ namespace Opm virtual void initPrimaryVariablesEvaluation() const override; - virtual void maybeDoGasLiftOptimization ( + virtual void gasLiftOptimizationStage1 ( WellState&, const Simulator&, - DeferredLogger& + DeferredLogger&, + GLiftProdWells &, + GLiftOptWells &, + GLiftWellStateMap & ) const override { // Not implemented yet } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index daba8d98f..61a2689c3 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -27,11 +27,13 @@ #include #endif -#include #include +#include +#include #include #include #include +#include #include #include @@ -46,6 +48,7 @@ #include #include +#include #include #include @@ -73,7 +76,10 @@ namespace Opm using typename Base::SparseMatrixAdapter; using typename Base::FluidState; using typename Base::RateVector; - using GasLiftHandler = Opm::GasLiftRuntime; + using typename Base::GasLiftSingleWell; + using typename Base::GLiftOptWells; + using typename Base::GLiftProdWells; + using typename Base::GLiftWellStateMap; using Base::numEq; using Base::numPhases; @@ -250,10 +256,13 @@ namespace Opm DeferredLogger& deferred_logger ) const; - virtual void maybeDoGasLiftOptimization ( + virtual void gasLiftOptimizationStage1 ( WellState& well_state, const Simulator& ebosSimulator, - DeferredLogger& deferred_logger + DeferredLogger& deferred_logger, + GLiftProdWells &prod_wells, + GLiftOptWells &glift_wells, + GLiftWellStateMap &state_map ) const override; bool checkGliftNewtonIterationIdxOk( @@ -399,6 +408,9 @@ namespace Opm // Enable GLIFT debug mode. This will enable output of logging messages. bool glift_debug = false; + // Optimize only wells under THP control + bool glift_optimize_only_thp_wells = true; + const EvalWell& getBhp() const; EvalWell getQs(const int comp_idx) const; diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 1b6e5d605..76d0d5043 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -275,7 +275,7 @@ namespace Opm return primary_variables_evaluation_[SFrac]; } } - else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { + else if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { if (FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && compIdx == Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx)) { return primary_variables_evaluation_[GFrac]; } @@ -300,7 +300,7 @@ namespace Opm well_fraction -= primary_variables_evaluation_[GFrac]; } - + return well_fraction; } @@ -1115,7 +1115,7 @@ namespace Opm if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { F[pu.phase_pos[Oil]] = 1.0; - + if (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) { F[pu.phase_pos[Water]] = primary_variables_[WFrac]; F[pu.phase_pos[Oil]] -= F[pu.phase_pos[Water]]; @@ -1213,7 +1213,7 @@ namespace Opm std::vector F(number_of_phases_, 0.0); double F_solvent = 0.0; if ( FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) ) { - const int oil_pos = pu.phase_pos[Oil]; + const int oil_pos = pu.phase_pos[Oil]; F[oil_pos] = 1.0; if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { @@ -1234,7 +1234,7 @@ namespace Opm } } else if ( FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx) ) { - const int water_pos = pu.phase_pos[Water]; + const int water_pos = pu.phase_pos[Water]; F[water_pos] = 1.0; if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { @@ -1244,7 +1244,7 @@ namespace Opm } } else if ( FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) ) { - const int gas_pos = pu.phase_pos[Gas]; + const int gas_pos = pu.phase_pos[Gas]; F[gas_pos] = 1.0; } @@ -2696,19 +2696,31 @@ namespace Opm gliftDebug("Optimization disabled in WellState", deferred_logger); return false; } - const int well_index = index_of_well_; - const Well::ProducerCMode& control_mode - = well_state.currentProductionControls()[well_index]; - if (control_mode != Well::ProducerCMode::THP ) { - gliftDebug("Not THP control", deferred_logger); + if (well_state.gliftCheckAlqOscillation(name())) { + gliftDebug("further optimization skipped due to oscillation in ALQ", + deferred_logger); return false; } + if (this->glift_optimize_only_thp_wells) { + const int well_index = index_of_well_; + const Well::ProducerCMode& control_mode + = well_state.currentProductionControls()[well_index]; + if (control_mode != Well::ProducerCMode::THP ) { + gliftDebug("Not THP control", deferred_logger); + return false; + } + } if (!checkGliftNewtonIterationIdxOk(ebos_simulator, deferred_logger)) { return false; } const int report_step_idx = ebos_simulator.episodeIndex(); const Opm::Schedule& schedule = ebos_simulator.vanguard().schedule(); const GasLiftOpt& glo = schedule.glo(report_step_idx); + if (!glo.has_well(name())) { + gliftDebug("Gas Lift not activated: WLIFTOPT is probably missing", + deferred_logger); + return false; + } auto increment = glo.gaslift_increment(); // NOTE: According to the manual: LIFTOPT, item 1, : // "Increment size for lift gas injection rate. Lift gas is @@ -2771,7 +2783,7 @@ namespace Opm { if (this->glift_debug) { const std::string message = fmt::format( - " GLIFT (DEBUG) : Well {} : {}", this->name(), msg); + " GLIFT (DEBUG) : SW : Well {} : {}", this->name(), msg); deferred_logger.info(message); } } @@ -2819,30 +2831,37 @@ namespace Opm ); } - template void StandardWell:: - maybeDoGasLiftOptimization( - WellState& well_state, - const Simulator& ebos_simulator, - Opm::DeferredLogger& deferred_logger) const + gasLiftOptimizationStage1( + WellState& well_state, + const Simulator& ebos_simulator, + Opm::DeferredLogger& deferred_logger, + GLiftProdWells &prod_wells, + GLiftOptWells &glift_wells, + GLiftWellStateMap &glift_state_map + //std::map &prod_wells + ) const { const auto& well = well_ecl_; if (well.isProducer()) { const auto& summary_state = ebos_simulator.vanguard().summaryState(); - const Well::ProducerCMode& current_control - = well_state.currentProductionControls()[this->index_of_well_]; - if ( this->Base::wellHasTHPConstraints(summary_state) - && current_control != Well::ProducerCMode::BHP ) { + if ( this->Base::wellHasTHPConstraints(summary_state) ) { if (doGasLiftOptimize(well_state, ebos_simulator, deferred_logger)) { - const auto& controls = well.productionControls(summary_state); - GasLiftHandler glift { - *this, ebos_simulator, summary_state, - deferred_logger, well_state, controls }; - glift.runOptimize(); + std::unique_ptr glift + = std::make_unique( + *this, ebos_simulator, summary_state, + deferred_logger, well_state); + auto state = glift->runOptimize(); + if (state) { + glift_state_map.insert({this->name(), std::move(state)}); + glift_wells.insert({this->name(), std::move(glift)}); + return; + } } } + prod_wells.insert({this->name(), this}); } } diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 2b10944f0..881e46902 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -41,6 +41,15 @@ #include #include #include +// NOTE: GasLiftSingleWell.hpp includes StandardWell.hpp which includes ourself +// (WellInterface.hpp), so we need to forward declare GasLiftSingleWell +// for it to be defined in this file. Similar for BlackoilWellModel +namespace Opm { + template class GasLiftSingleWell; + template class BlackoilWellModel; +} +#include +#include #include #include @@ -83,6 +92,11 @@ namespace Opm using MaterialLaw = GetPropType; using SparseMatrixAdapter = GetPropType; using RateVector = GetPropType; + using GasLiftSingleWell = Opm::GasLiftSingleWell; + using GLiftOptWells = typename Opm::BlackoilWellModel::GLiftOptWells; + using GLiftProdWells = typename Opm::BlackoilWellModel::GLiftProdWells; + using GLiftWellStateMap = + typename Opm::BlackoilWellModel::GLiftWellStateMap; static const int numEq = Indices::numEq; static const int numPhases = Indices::numPhases; @@ -174,10 +188,13 @@ namespace Opm Opm::DeferredLogger& deferred_logger ) = 0; - virtual void maybeDoGasLiftOptimization ( + virtual void gasLiftOptimizationStage1 ( WellState& well_state, const Simulator& ebosSimulator, - DeferredLogger& deferred_logger + DeferredLogger& deferred_logger, + GLiftProdWells& prod_wells, + GLiftOptWells& glift_wells, + GLiftWellStateMap& state_map ) const = 0; void updateWellTestState(const WellState& well_state, @@ -324,6 +341,7 @@ namespace Opm WellState& well_state, Opm::DeferredLogger& deferred_logger); + const PhaseUsage& phaseUsage() const; protected: @@ -435,8 +453,6 @@ namespace Opm bool changed_to_stopped_this_step_ = false; - const PhaseUsage& phaseUsage() const; - int flowPhaseToEbosCompIdx( const int phaseIdx ) const; int flowPhaseToEbosPhaseIdx( const int phaseIdx ) const; diff --git a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp index 137451278..16e8d71b2 100644 --- a/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp +++ b/opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp @@ -1283,10 +1283,66 @@ namespace Opm this->current_alq_[name] = value; } + bool gliftCheckAlqOscillation(const std::string &name) const { + if ((this->alq_increase_count_.count(name) == 1) && + (this->alq_decrease_count_.count(name) == 1)) + { + if ((this->alq_increase_count_.at(name) >= 1) && + (this->alq_decrease_count_.at(name) >= 1)) + { + return true; + } + } + return false; + } + + int gliftGetAlqDecreaseCount(const std::string &name) { + if (this->alq_decrease_count_.count(name) == 0) { + return 0; + } + else { + return this->alq_decrease_count_[name]; + } + } + + int gliftGetAlqIncreaseCount(const std::string &name) { + if (this->alq_increase_count_.count(name) == 0) { + return 0; + } + else { + return this->alq_increase_count_[name]; + } + } + + void gliftUpdateAlqIncreaseCount(const std::string &name, bool increase) { + if (increase) { + if (this->alq_increase_count_.count(name) == 0) { + this->alq_increase_count_[name] = 1; + } + else { + this->alq_increase_count_[name]++; + } + } + else { + if (this->alq_decrease_count_.count(name) == 0) { + this->alq_decrease_count_[name] = 1; + } + else { + this->alq_decrease_count_[name]++; + } + } + } + bool gliftOptimizationEnabled() const { return do_glift_optimization_; } + void gliftTimeStepInit() { + this->alq_increase_count_.clear(); + this->alq_decrease_count_.clear(); + disableGliftOptimization(); + } + void disableGliftOptimization() { do_glift_optimization_ = false; } @@ -1325,6 +1381,8 @@ namespace Opm std::map group_grat_target_from_sales; std::map current_alq_; std::map default_alq_; + std::map alq_increase_count_; + std::map alq_decrease_count_; bool do_glift_optimization_; std::vector perfRateSolvent_;