From 434640fdf5ab4ae507aa5c445dd7c906f47449f3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?H=C3=A5kon=20H=C3=A6gland?= Date: Wed, 3 Mar 2021 13:59:26 +0100 Subject: [PATCH] Implements gas lift optimization for groups. Extends PR #2824 to include support for GLIFTOPT (item 2, maximum lift gas supply for a group) and group production constraints. The optimization is split into two phases. First the wells are optimized separately (as in PR #2824). In this phase LIFTOPT and WLIFTOPT constraints (e.g. maxmimum lift gas injection for a well, minimum economic gradient) are considered together with well production constraints. Then, in the next phase the wells are optimized in groups. Here, the ALQ distribution from the first phase is used as a starting point. If a group has any production rate constraints, and/or a limit on its total rate of lift gas supply, lift gas is redistributed to the wells that gain the most benefit from it by considering which wells that currently has the largest weighted incremental gradient (i.e. increase in oil rate compared to increase in ALQ). --- opm/simulators/flow/BlackoilModelEbos.hpp | 4 +- opm/simulators/wells/BlackoilWellModel.hpp | 24 +- .../wells/BlackoilWellModel_impl.hpp | 59 +- opm/simulators/wells/GasLiftRuntime.hpp | 153 -- opm/simulators/wells/GasLiftRuntime_impl.hpp | 816 ---------- opm/simulators/wells/GasLiftSingleWell.hpp | 235 +++ .../wells/GasLiftSingleWell_impl.hpp | 1370 +++++++++++++++++ opm/simulators/wells/GasLiftStage2.hpp | 225 +++ opm/simulators/wells/GasLiftStage2_impl.hpp | 1027 ++++++++++++ opm/simulators/wells/GasLiftWellState.hpp | 82 + .../wells/GasLiftWellState_impl.hpp | 18 + opm/simulators/wells/MultisegmentWell.hpp | 11 +- opm/simulators/wells/StandardWell.hpp | 20 +- opm/simulators/wells/StandardWell_impl.hpp | 71 +- opm/simulators/wells/WellInterface.hpp | 24 +- .../wells/WellStateFullyImplicitBlackoil.hpp | 58 + 16 files changed, 3185 insertions(+), 1012 deletions(-) delete mode 100644 opm/simulators/wells/GasLiftRuntime.hpp delete mode 100644 opm/simulators/wells/GasLiftRuntime_impl.hpp create mode 100644 opm/simulators/wells/GasLiftSingleWell.hpp create mode 100644 opm/simulators/wells/GasLiftSingleWell_impl.hpp create mode 100644 opm/simulators/wells/GasLiftStage2.hpp create mode 100644 opm/simulators/wells/GasLiftStage2_impl.hpp create mode 100644 opm/simulators/wells/GasLiftWellState.hpp create mode 100644 opm/simulators/wells/GasLiftWellState_impl.hpp 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_;