diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 2658e1f55..f76d6406f 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -60,6 +60,7 @@ list (APPEND MAIN_SOURCE_FILES opm/simulators/utils/ParallelRestart.cpp opm/simulators/wells/ALQState.cpp opm/simulators/wells/BlackoilWellModelGeneric.cpp + opm/simulators/wells/GasLiftGroupInfo.cpp opm/simulators/wells/GasLiftSingleWellGeneric.cpp opm/simulators/wells/GasLiftStage2.cpp opm/simulators/wells/GlobalWellInfo.cpp diff --git a/opm/simulators/wells/BlackoilWellModel.hpp b/opm/simulators/wells/BlackoilWellModel.hpp index 4929f8224..04208ddea 100644 --- a/opm/simulators/wells/BlackoilWellModel.hpp +++ b/opm/simulators/wells/BlackoilWellModel.hpp @@ -31,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -50,6 +51,9 @@ #include #include #include +#include +#include +#include #include #include #include @@ -100,6 +104,12 @@ namespace Opm { using RateVector = GetPropType; using GlobalEqVector = GetPropType; using SparseMatrixAdapter = GetPropType; + using GLiftOptWells = typename BlackoilWellModelGeneric::GLiftOptWells; + using GLiftProdWells = typename BlackoilWellModelGeneric::GLiftProdWells; + using GLiftWellStateMap = + typename BlackoilWellModelGeneric::GLiftWellStateMap; + using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells; + using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; typedef typename BaseAuxiliaryModule::NeighborSet NeighborSet; @@ -266,6 +276,7 @@ namespace Opm { void initPrimaryVariablesEvaluation() const; void updateWellControls(DeferredLogger& deferred_logger, const bool checkGroupControls); WellInterfacePtr getWell(const std::string& well_name) const; + void initGliftEclWellMap(GLiftEclWells &ecl_well_map); protected: Simulator& ebosSimulator_; @@ -360,6 +371,12 @@ namespace Opm { void maybeDoGasLiftOptimize(DeferredLogger& deferred_logger); + bool checkDoGasLiftOptimization(DeferredLogger& deferred_logger); + + void gasLiftOptimizationStage1(DeferredLogger& deferred_logger, + GLiftProdWells &prod_wells, GLiftOptWells &glift_wells, + GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map); + void extractLegacyCellPvtRegionIndex_(); void extractLegacyDepth_(); diff --git a/opm/simulators/wells/BlackoilWellModel_impl.hpp b/opm/simulators/wells/BlackoilWellModel_impl.hpp index c30487da5..eecafd8b0 100644 --- a/opm/simulators/wells/BlackoilWellModel_impl.hpp +++ b/opm/simulators/wells/BlackoilWellModel_impl.hpp @@ -851,19 +851,165 @@ namespace Opm { BlackoilWellModel:: maybeDoGasLiftOptimize(DeferredLogger& deferred_logger) { - this->wellState().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( - this->wellState(), ebosSimulator_, deferred_logger, - prod_wells, glift_wells, state_map); + if (checkDoGasLiftOptimization(deferred_logger)) { + GLiftOptWells glift_wells; + GLiftProdWells prod_wells; + GLiftWellStateMap state_map; + // NOTE: To make GasLiftGroupInfo (see below) independent of the TypeTag + // associated with *this (i.e. BlackoilWellModel) we observe + // that GasLiftGroupInfo's only dependence on *this is that it needs to + // access the eclipse Wells in the well container (the eclipse Wells + // themselves are independent of the TypeTag). + // Hence, we extract them from the well container such that we can pass + // them to the GasLiftGroupInfo constructor. + GLiftEclWells ecl_well_map; + initGliftEclWellMap(ecl_well_map); + GasLiftGroupInfo group_info { + ecl_well_map, + ebosSimulator_.vanguard().schedule(), + ebosSimulator_.vanguard().summaryState(), + ebosSimulator_.episodeIndex(), + ebosSimulator_.model().newtonMethod().numIterations(), + phase_usage_, + deferred_logger, + this->wellState() + }; + group_info.initialize(ebosSimulator_.vanguard().grid().comm()); + gasLiftOptimizationStage1( + deferred_logger, prod_wells, glift_wells, group_info, state_map); + gasLiftOptimizationStage2( + deferred_logger, prod_wells, glift_wells, state_map, + ebosSimulator_.episodeIndex()); + if (this->glift_debug) gliftDebugShowALQ(deferred_logger); } - gasLiftOptimizationStage2(deferred_logger, prod_wells, glift_wells, state_map, ebosSimulator_.episodeIndex()); - if (this->glift_debug) gliftDebugShowALQ(deferred_logger); - this->wellState().disableGliftOptimization(); + } + + template + void + BlackoilWellModel:: + gasLiftOptimizationStage1(DeferredLogger& deferred_logger, + GLiftProdWells &prod_wells, GLiftOptWells &glift_wells, + GasLiftGroupInfo &group_info, GLiftWellStateMap &state_map) + { + auto comm = ebosSimulator_.vanguard().grid().comm(); + int num_procs = comm.size(); + // NOTE: Gas lift optimization stage 1 seems to be difficult + // to do in parallel since the wells are optimized on different + // processes and each process needs to know the current ALQ allocated + // to each group it is a memeber of in order to check group limits and avoid + // allocating more ALQ than necessary. (Surplus ALQ is removed in + // stage 2). In stage1, as each well is adding ALQ, the current group ALQ needs + // to be communicated to the other processes. But there is no common + // synchronization point that all process will reach in the + // runOptimizeLoop_() in GasLiftSingleWell.cpp. + // + // TODO: Maybe a better solution could be invented by distributing + // wells according to certain parent groups. Then updated group rates + // might not have to be communicated to the other processors. + + // Currently, the best option seems to be to run this part sequentially + // (not in parallel). + // + // TODO: The simplest approach seems to be if a) one process could take + // ownership of all the wells (the union of all the wells in the + // well_container_ of each process) then this process could do the + // optimization, while the other processes could wait for it to + // finish (e.g. comm.barrier()), or alternatively, b) if all + // processes could take ownership of all the wells. Then there + // would be no need for synchronization here.. + // + for (int i = 0; i< num_procs; i++) { + int num_rates_to_sync = 0; // communication variable + GLiftSyncGroups groups_to_sync; + if (comm.rank() == i) { + // Run stage1: Optimize single wells while also checking group limits + for (const auto& well : well_container_) { + // NOTE: Only the wells in "group_info" needs to be optimized + if (group_info.hasWell(well->name())) { + well->gasLiftOptimizationStage1( + this->wellState(), ebosSimulator_, deferred_logger, + prod_wells, glift_wells, state_map, + group_info, groups_to_sync); + } + } + num_rates_to_sync = groups_to_sync.size(); + } + // Since "group_info" is not used in stage2, there is no need to + // communicate rates if this is the last iteration... + if (i == (num_procs - 1)) + break; + num_rates_to_sync = comm.sum(num_rates_to_sync); + if (num_rates_to_sync > 0) { + std::vector group_indexes; + group_indexes.reserve(num_rates_to_sync); + std::vector group_alq_rates; + group_alq_rates.reserve(num_rates_to_sync); + std::vector group_oil_rates; + group_oil_rates.reserve(num_rates_to_sync); + std::vector group_gas_rates; + group_gas_rates.reserve(num_rates_to_sync); + if (comm.rank() == i) { + for (auto idx : groups_to_sync) { + auto [oil_rate, gas_rate, alq] = group_info.getRates(idx); + group_indexes.push_back(idx); + group_oil_rates.push_back(oil_rate); + group_gas_rates.push_back(gas_rate); + group_alq_rates.push_back(alq); + } + } + // TODO: We only need to broadcast to processors with index + // j > i since we do not use the "group_info" in stage 2. In + // this case we should use comm.send() and comm.receive() + // instead of comm.broadcast() to communicate with specific + // processes, and these processes only need to receive the + // data if they are going to check the group rates in stage1 + // Another similar idea is to only communicate the rates to + // process j = i + 1 + comm.broadcast(group_indexes.data(), num_rates_to_sync, i); + comm.broadcast(group_oil_rates.data(), num_rates_to_sync, i); + comm.broadcast(group_gas_rates.data(), num_rates_to_sync, i); + comm.broadcast(group_alq_rates.data(), num_rates_to_sync, i); + if (comm.rank() != i) { + for (int j=0; j + void + BlackoilWellModel:: + initGliftEclWellMap(GLiftEclWells &ecl_well_map) + { + for ( const auto& well: well_container_ ) { + ecl_well_map.try_emplace( + well->name(), &(well->wellEcl()), well->indexOfWell()); + } + } + + template + bool + BlackoilWellModel:: + checkDoGasLiftOptimization(Opm::DeferredLogger& deferred_logger) + { + gliftDebug("checking if GLIFT should be done..", deferred_logger); + /* + std::size_t num_procs = ebosSimulator_.gridView().comm().size(); + if (num_procs > 1u) { + const std::string msg = fmt::format(" GLIFT: skipping optimization. " + "Parallel run not supported yet: num procs = {}", num_procs); + deferred_logger.warning(msg); + return false; + } + */ + if (!(this->wellState().gliftOptimizationEnabled())) { + gliftDebug("Optimization disabled in WellState", deferred_logger); + return false; + } + return true; } template diff --git a/opm/simulators/wells/GasLiftGroupInfo.cpp b/opm/simulators/wells/GasLiftGroupInfo.cpp new file mode 100644 index 000000000..ca1df132a --- /dev/null +++ b/opm/simulators/wells/GasLiftGroupInfo.cpp @@ -0,0 +1,373 @@ +/* + 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 + +namespace Opm { + +GasLiftGroupInfo:: +GasLiftGroupInfo( + GLiftEclWells &ecl_wells, + const Schedule &schedule, + const SummaryState &summary_state, + const int report_step_idx, + const int iteration_idx, + const PhaseUsage &phase_usage, + DeferredLogger &deferred_logger, + WellState &well_state +) : + ecl_wells_{ecl_wells}, + schedule_{schedule}, + summary_state_{summary_state}, + report_step_idx_{report_step_idx}, + iteration_idx_{iteration_idx}, + phase_usage_{phase_usage}, + deferred_logger_{deferred_logger}, + well_state_{well_state}, + glo_{schedule_.glo(report_step_idx_)}, + debug{false} +{ + +} + +/**************************************** + * Public methods in alphabetical order + ****************************************/ + +double +GasLiftGroupInfo:: +alqRate(const std::string& group_name) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + return group_rate.alq(); +} + +int +GasLiftGroupInfo:: +getGroupIdx(const std::string& group_name) +{ + return this->group_idx_.at(group_name); +} + +double +GasLiftGroupInfo:: +gasRate(const std::string& group_name) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + return group_rate.gasRate(); +} + +std::optional +GasLiftGroupInfo:: +gasTarget(const std::string& group_name) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + return group_rate.gasTarget(); +} + +std::tuple +GasLiftGroupInfo:: +getRates(int group_idx) +{ + const auto& group_name = groupIdxToName(group_idx); + auto& rates = this->group_rate_map_.at(group_name); + return std::make_tuple(rates.oilRate(), rates.gasRate(), rates.alq()); +} + +std::vector>& +GasLiftGroupInfo:: +getWellGroups(const std::string& well_name) +{ + assert(this->well_group_map_.count(well_name) == 1); + return this->well_group_map_[well_name]; +} + +const std::string& +GasLiftGroupInfo:: +groupIdxToName(int group_idx) +{ + const std::string *group_name = nullptr; + // TODO: An alternative to the below loop is to set up a reverse map from idx -> + // string, then we could in theory do faster lookup here.. + for (const auto& [key, value] : this->group_idx_) { + if (value == group_idx) { + // NOTE: it is assumed that the mapping from name->idx is one-to-one + // so there can only be one idx with a given group name. + group_name = &key; + break; + } + } + // the caller is responsible for providing a valid idx, so group_name + // cannot be nullptr here.. + assert(group_name); + return *group_name; +} + +bool +GasLiftGroupInfo:: +hasWell(const std::string& well_name) +{ + return this->well_group_map_.count(well_name) == 1; +} + + +std::optional +GasLiftGroupInfo:: +maxAlq(const std::string& group_name) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + return group_rate.maxAlq(); +} + +double +GasLiftGroupInfo:: +oilRate(const std::string &group_name) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + return group_rate.oilRate(); +} + +std::optional +GasLiftGroupInfo:: +oilTarget(const std::string &group_name) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + return group_rate.oilTarget(); +} + +void +GasLiftGroupInfo:: +update( + const std::string &group_name, double delta_oil, double delta_gas, double delta_alq) +{ + auto& group_rate = this->group_rate_map_.at(group_name); + group_rate.update(delta_oil, delta_gas, delta_alq); +} + +void +GasLiftGroupInfo:: +updateRate(int idx, double oil_rate, double gas_rate, double alq) +{ + const auto& group_name = groupIdxToName(idx); + auto& rates = this->group_rate_map_.at(group_name); + rates.assign(oil_rate, gas_rate, alq); +} + +/**************************************** + * Private methods in alphabetical order + ****************************************/ + + +bool +GasLiftGroupInfo:: +checkDoGasLiftOptimization_(const std::string &well_name) +{ + if (this->well_state_.gliftCheckAlqOscillation(well_name)) { + displayDebugMessage_( + "further optimization skipped due to oscillation in ALQ", well_name); + return false; + } + if (this->optimize_only_thp_wells_) { + auto itr = this->ecl_wells_.find(well_name); + if (itr != this->ecl_wells_.end()) { + //const Well *well = (itr->second).first; + //assert(well); // Should never be nullptr + const int index = (itr->second).second; + const Well::ProducerCMode& control_mode + = this->well_state_.currentProductionControl(index); + if (control_mode != Well::ProducerCMode::THP ) { + displayDebugMessage_("Not THP control. Skipping.", well_name); + return false; + } + } + else { + // well_name is not present in the well_model's well container + return false; + } + } + if (!checkNewtonIterationIdxOk_(well_name)) { + return false; + } + if (!this->glo_.has_well(well_name)) { + displayDebugMessage_( + "Gas Lift not activated: WLIFTOPT is probably missing", well_name); + return false; + } + auto increment = this->glo_.gaslift_increment(); + // NOTE: According to the manual: 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." + if (increment <= 0) { + if (this->debug) { + const std::string msg = fmt::format( + "Gas Lift switched off in LIFTOPT item 1 due to non-positive " + "value: {}", increment); + displayDebugMessage_(msg, well_name); + } + return false; + } + else { + return true; + } +} + +bool +GasLiftGroupInfo:: +checkNewtonIterationIdxOk_(const std::string &well_name) +{ + if (this->glo_.all_newton()) { + const int nupcol = this->schedule_[this->report_step_idx_].nupcol(); + if (this->debug) { + const std::string msg = fmt::format( + "LIFTOPT item4 == YES, it = {}, nupcol = {} --> GLIFT optimize = {}", + this->iteration_idx_, + nupcol, + ((this->iteration_idx_ <= nupcol) ? "TRUE" : "FALSE")); + displayDebugMessage_(msg, well_name); + } + return this->iteration_idx_ <= nupcol; + } + else { + if (this->debug) { + const std::string msg = fmt::format( + "LIFTOPT item4 == NO, it = {} --> GLIFT optimize = {}", + this->iteration_idx_, + ((this->iteration_idx_ == 1) ? "TRUE" : "FALSE")); + displayDebugMessage_(msg, well_name); + } + return this->iteration_idx_ == 1; + } +} + +void +GasLiftGroupInfo:: +displayDebugMessage_(const std::string &msg) +{ + if (this->debug) { + const std::string message = fmt::format( + " GLIFT (DEBUG) : Init group info : {}", msg); + this->deferred_logger_.info(message); + } +} + +void +GasLiftGroupInfo:: +displayDebugMessage_(const std::string &msg, const std::string &well_name) +{ + if (this->debug) { + const std::string message = fmt::format( + " GLIFT (DEBUG) : Init group info : Well {} : {}", + well_name, msg); + this->deferred_logger_.info(message); + } +} + + +std::pair +GasLiftGroupInfo:: +getProducerWellRates_(int well_index) +{ + const auto& pu = this->phase_usage_; + auto oil_rate = + -this->well_state_.wellRates(well_index)[pu.phase_pos[Oil]]; + auto gas_rate = + -this->well_state_.wellRates(well_index)[pu.phase_pos[Gas]]; + return {oil_rate, gas_rate}; +} + +void +GasLiftGroupInfo:: +initializeWell2GroupMapRecursive_( + const Group &group, + std::vector &group_names, + std::vector &group_efficiency, + double cur_efficiency) +{ + double gfac = group.getGroupEfficiencyFactor(); + cur_efficiency = gfac * cur_efficiency; + for (auto &item : group_efficiency) { + item *= gfac; + } + if (this->group_rate_map_.count(group.name()) == 1) { + // extract the subset of groups that has limits or targets that can affect + // gas lift optimization. + group_names.push_back(group.name()); + group_efficiency.push_back(gfac); + } + if (group.wellgroup()) { + for (const std::string& well_name : group.wells()) { + // TODO: can the same well be memember of two different groups + // (on the same recursion level) ? + assert(this->well_group_map_.count(well_name) == 0); + if (checkDoGasLiftOptimization_(well_name)) { + const auto &well = this->schedule_.getWell( + well_name, this->report_step_idx_); + double wfac = well.getEfficiencyFactor(); + auto [itr, success] = this->well_group_map_.insert( + {well_name, /*empty vector*/ {}}); + assert(success); + auto &vec = itr->second; + assert(group_names.size() == group_efficiency.size()); + auto iter2 = group_efficiency.begin(); + for (auto iter1 = group_names.begin(); + iter1 != group_names.end(); ++iter1) + { + double efficiency = (*iter2) * wfac; + vec.emplace_back(/*group_name=*/*iter1, efficiency); + ++iter2; + } + } + } + } + else { + 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_); + initializeWell2GroupMapRecursive_( + sub_group, group_names, group_efficiency, cur_efficiency); + } + } + if (this->group_rate_map_.count(group.name()) == 1) { + group_names.pop_back(); + group_efficiency.pop_back(); + } +} + + + +// TODO: It would be more efficient if the group idx map was build once +// per time step (or better: once per report step) and saved e.g. in +// the well state object, instead of rebuilding here for each of +// NUPCOL well iteration for each time step. +void +GasLiftGroupInfo:: +updateGroupIdxMap_(const std::string &group_name) +{ + if (this->group_idx_.count(group_name) == 0) { + //auto [itr, success] = + this->group_idx_.try_emplace(group_name, this->next_group_idx_); + this->next_group_idx_++; + } +} + +} // namespace Opm diff --git a/opm/simulators/wells/GasLiftGroupInfo.hpp b/opm/simulators/wells/GasLiftGroupInfo.hpp new file mode 100644 index 000000000..f478684c2 --- /dev/null +++ b/opm/simulators/wells/GasLiftGroupInfo.hpp @@ -0,0 +1,274 @@ +/* + 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_GROUP_INFO_HEADER_INCLUDED +#define OPM_GASLIFT_GROUP_INFO_HEADER_INCLUDED + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +namespace Opm +{ +class GasLiftGroupInfo +{ + class GroupRates; + // NOTE: In the Well2GroupMap below, in the std::vector value we store + // pairs of group names and efficiency factors. The efficiency factors + // are the product of the wells efficiency factor and all the efficiency + // factors of the child groups of the group all the way down + // to the well group. + using Well2GroupMap = + std::map>>; + using GroupRateMap = + std::map; + using GroupIdxMap = std::map; + + // TODO: same definition with WellInterface, and + // WellState 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; +public: + using GLiftEclWells = std::map>; + GasLiftGroupInfo( + GLiftEclWells& ecl_wells, + const Schedule& schedule, + const SummaryState& summary_state, + const int report_step_idx, + const int iteration_idx, + const PhaseUsage& phase_usage, + DeferredLogger& deferred_logger, + WellState& well_state); + std::vector>& getWellGroups( + const std::string& well_name); + + // TODO: See comment below for initializeGroupRatesRecursive_() for why + // the implementation of initialize() is kept here in the header file instead + // of in the .cpp file... + template + void + initialize(const Comm& comm) + { + const auto& group = this->schedule_.getGroup("FIELD", this->report_step_idx_); + initializeGroupRatesRecursive_(comm, group); + std::vector group_names; + std::vector group_efficiency; + initializeWell2GroupMapRecursive_( + group, group_names, group_efficiency, /*current efficiency=*/1.0); + } + + double alqRate(const std::string& group_name); + double gasRate(const std::string& group_name); + int getGroupIdx(const std::string& group_name); + std::tuple getRates(int group_idx); + std::optional gasTarget(const std::string& group_name); + const std::string& groupIdxToName(int group_idx); + bool hasWell(const std::string& well_name); + std::optional maxAlq(const std::string& group_name); + double oilRate(const std::string& group_name); + std::optional oilTarget(const std::string& group_name); + void update(const std::string& well_name, + double delta_oil, double delta_gas, double delta_alq); + void updateRate(int idx, double oil_rate, double gas_rate, double alq); + const Well2GroupMap& wellGroupMap() { return well_group_map_; } +private: + bool checkDoGasLiftOptimization_(const std::string& well_name); + bool checkNewtonIterationIdxOk_(const std::string& well_name); + void displayDebugMessage_(const std::string& msg); + void displayDebugMessage_(const std::string& msg, const std::string& well_name); + std::pair getProducerWellRates_(const int index); + void initializeWell2GroupMapRecursive_( + const Group& group, std::vector& group_names, + std::vector& group_efficiency, double cur_efficiency); + void updateGroupIdxMap_(const std::string& group_name); + + // TODO: I first tried to pass the MPI Communicator as a constructor argument + // to this class (GasLiftGroupInfo) similar to what is done for + // GasLiftStage2 (see GasLiftStage2.hpp), hower this did not work for this + // class since we are also constructing a GasLiftGroupInfo object in the + // test file test1_glift.cpp and when the linker tries to find a definition + // of the GasLiftGroupInfo(...) constructor in libopmsimulators.a, + // the template type of the MPI communicator (Dune::Communication<..>) + // is not of the same type as the one needed by the test case. + // The test case needs Dune::Communication, whereas + // the one in libopmsimulators.a is Dune::Communication. + // + // The error I got from the linker is: + // + // /bin/ld: CMakeFiles/test_glift1.dir/tests/test_glift1.cpp.o: + // in function `G1::test_method()': + // test_glift1.cpp:(.text+0x15b36): undefined reference to + // `Opm::GasLiftGroupInfo::GasLiftGroupInfo(....) + // + // to work around this issue this function is templetized in terms of Comm + // here in the header file instead of having it in the .cpp file. + // (thanks to Tor Harald S. for the suggestion) + template + std::tuple + initializeGroupRatesRecursive_(const Comm& comm, const Group &group) + { + double oil_rate = 0.0; + double gas_rate = 0.0; + double alq = 0.0; + if (group.wellgroup()) { + for (const std::string& well_name : group.wells()) { + // NOTE: we cannot simply use: + // + // const auto &well = + // this->schedule_.getWell(well_name, this->report_step_idx_); + // + // since the well may not be active (present in the well container) + auto itr = this->ecl_wells_.find(well_name); + if (itr != this->ecl_wells_.end()) { + const Well *well = (itr->second).first; + assert(well); // Should never be nullptr + const int index = (itr->second).second; + if (well->isProducer()) { + auto [sw_oil_rate, sw_gas_rate] = getProducerWellRates_(index); + auto sw_alq = this->well_state_.getALQ(well_name); + double factor = well->getEfficiencyFactor(); + oil_rate += (factor * sw_oil_rate); + gas_rate += (factor * sw_gas_rate); + alq += (factor * sw_alq); + } + } + } + // These sums needs to be communictated + oil_rate = comm.sum(oil_rate); + gas_rate = comm.sum(gas_rate); + alq = comm.sum(alq); + } + else { + 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_); + auto [sg_oil_rate, sg_gas_rate, sg_alq] + = initializeGroupRatesRecursive_(comm, sub_group); + const auto gefac = sub_group.getGroupEfficiencyFactor(); + oil_rate += (gefac * sg_oil_rate); + gas_rate += (gefac * sg_gas_rate); + alq += (gefac * sg_alq); + } + } + std::optional oil_target, gas_target, max_total_gas, max_alq; + const auto controls = group.productionControls(this->summary_state_); + if (group.has_control(Group::ProductionCMode::ORAT)) { + oil_target = controls.oil_target; + } + if (group.has_control(Group::ProductionCMode::GRAT)) { + gas_target = controls.gas_target; + } + if (this->glo_.has_group(group.name())) { + const auto &gl_group = this->glo_.group(group.name()); + max_alq = gl_group.max_lift_gas(); + max_total_gas = gl_group.max_total_gas(); + } + if (oil_target || gas_target || max_total_gas || max_alq) { + updateGroupIdxMap_(group.name()); + this->group_rate_map_.try_emplace(group.name(), + oil_rate, gas_rate, alq, oil_target, gas_target, max_total_gas, max_alq); + } + return std::make_tuple(oil_rate, gas_rate, alq); + } + + class GroupRates { + public: + GroupRates( double oil_rate, double gas_rate, double alq, + std::optional oil_target, + std::optional gas_target, + std::optional total_gas, + std::optional max_alq + ) : + oil_rate_{oil_rate}, + gas_rate_{gas_rate}, + alq_{alq}, + oil_target_{oil_target}, + gas_target_{gas_target}, + total_gas_{total_gas}, + max_alq_{max_alq} + {} + double alq() const { return alq_; } + void assign(double oil_rate, double gas_rate, double alq) + { + oil_rate_ = oil_rate; + gas_rate_ = gas_rate; + alq_ = alq; + } + double gasRate() const { return gas_rate_; } + std::optional gasTarget() const { return gas_target_; } + std::optional maxAlq() const { return max_alq_; } + double oilRate() const { return oil_rate_; } + std::optional oilTarget() const { return oil_target_; } + void update(double delta_oil, double delta_gas, double delta_alq) + { + oil_rate_ += delta_oil; + gas_rate_ += delta_gas; + alq_ += delta_alq; + } + private: + double oil_rate_; + double gas_rate_; + double alq_; + std::optional oil_target_; + std::optional gas_target_; + std::optional total_gas_; + std::optional max_alq_; + }; + + GLiftEclWells &ecl_wells_; + const Schedule &schedule_; + const SummaryState &summary_state_; + const int report_step_idx_; + const int iteration_idx_; + const PhaseUsage &phase_usage_; + DeferredLogger &deferred_logger_; + WellState &well_state_; + const GasLiftOpt& glo_; + GroupRateMap group_rate_map_; + Well2GroupMap well_group_map_; + bool debug; + GroupIdxMap group_idx_; + int next_group_idx_ = 0; + // Optimize only wells under THP control + bool optimize_only_thp_wells_ = true; + +}; + +} // namespace Opm + +#endif // OPM_GASLIFT_GROUP_INFO_INCLUDED diff --git a/opm/simulators/wells/GasLiftSingleWell.hpp b/opm/simulators/wells/GasLiftSingleWell.hpp index 5a744dde2..e08055daa 100644 --- a/opm/simulators/wells/GasLiftSingleWell.hpp +++ b/opm/simulators/wells/GasLiftSingleWell.hpp @@ -31,6 +31,7 @@ namespace Opm { } #include #include +#include #include #include @@ -44,6 +45,7 @@ namespace Opm { using Simulator = GetPropType; using StdWell = StandardWell; + using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; public: GasLiftSingleWell( @@ -51,7 +53,9 @@ namespace Opm const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, - WellState &well_state + WellState &well_state, + GasLiftGroupInfo &group_info, + GLiftSyncGroups &sync_groups ); const WellInterfaceGeneric &getStdWell() const override { return std_well_; } diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp index de0fb41b2..94b11adcd 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.cpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.cpp @@ -35,17 +35,23 @@ namespace Opm { -GasLiftSingleWellGeneric::GasLiftSingleWellGeneric(DeferredLogger& deferred_logger, - WellState& well_state, - const Well& ecl_well, - const SummaryState& summary_state, - const Schedule& schedule, - const int report_step_idx) - : deferred_logger_(deferred_logger) - , well_state_(well_state) - , ecl_well_(ecl_well) - , summary_state_(summary_state) - , controls_(ecl_well_.productionControls(summary_state_)) +GasLiftSingleWellGeneric::GasLiftSingleWellGeneric( + DeferredLogger& deferred_logger, + WellState& well_state, + const Well& ecl_well, + const SummaryState& summary_state, + GasLiftGroupInfo &group_info, + const Schedule& schedule, + const int report_step_idx, + GLiftSyncGroups &sync_groups +) : + deferred_logger_{deferred_logger} + , well_state_{well_state} + , ecl_well_{ecl_well} + , summary_state_{summary_state} + , group_info_{group_info} + , sync_groups_{sync_groups} + , controls_{ecl_well_.productionControls(summary_state_)} , num_phases_{well_state_.numPhases()} , debug_{false} // extra debugging output , debug_limit_increase_decrease_{false} @@ -106,33 +112,28 @@ calcIncOrDecGradient(double oil_rate, double gas_rate, double alq, bool increase } } -bool +std::unique_ptr GasLiftSingleWellGeneric:: -computeInitialWellRates_(std::vector& potentials) +runOptimize(const int iteration_idx) { - 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); + std::unique_ptr state; + if (this->optimize_) { + if (this->debug_limit_increase_decrease_) { + state = runOptimize1_(); } - 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); + else { + state = runOptimize2_(); + } + if (state) { + if (state->increase()) { + double alq = state->alq(); + if (this->debug_) + logSuccess_(alq, iteration_idx); + this->well_state_.setALQ(this->well_name_, alq); + } } - return true; - } - else { - displayDebugMessage_("Aborting optimization."); - return false; } + return state; } /**************************************** @@ -229,9 +230,10 @@ checkInitialALQmodified_(double alq, double initial_alq) const bool GasLiftSingleWellGeneric:: -checkWellRatesViolated_(std::vector& potentials, - const std::function& callback, - bool increase) +checkWellRatesViolated_( + std::vector& potentials, + const std::function& callback, + bool increase) { if (!increase) { auto oil_rate = -potentials[this->oil_pos_]; @@ -276,6 +278,35 @@ checkWellRatesViolated_(std::vector& potentials, return false; } +bool +GasLiftSingleWellGeneric:: +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; + } +} + void GasLiftSingleWellGeneric:: debugCheckNegativeGradient_(double grad, double alq, double new_alq, double oil_rate, @@ -566,6 +597,365 @@ logSuccess_(double alq, const int iteration_idx) this->deferred_logger_.info(message); } +std::tuple +GasLiftSingleWellGeneric:: +maybeAdjustALQbeforeOptimizeLoop_( + bool increase, double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, + std::vector &potentials) +{ + double orig_alq = alq; + if (this->debug_) { + const std::string msg = fmt::format("initial ALQ: {}", alq); + displayDebugMessage_(msg); + } + 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, alq) = + reduceALQtoOilTarget_( + 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, alq) = + increaseALQtoPositiveOilRate_(alq, oil_rate, + gas_rate, oil_is_limited, gas_is_limited, potentials); + } + if (increase && (this->min_alq_> 0) && (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, alq) = + increaseALQtoMinALQ_(alq, oil_rate, gas_rate, + oil_is_limited, gas_is_limited, potentials); + } + } + if (this->debug_ && (orig_alq != alq)) { + const std::string msg = fmt::format("adjusted ALQ to: {}", alq); + displayDebugMessage_(msg); + } + return std::make_tuple(oil_rate, gas_rate, alq, oil_is_limited, gas_is_limited); +} + +std::tuple +GasLiftSingleWellGeneric:: +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_) { + assert( alq <= orig_alq ); + if (alq < orig_alq) { + // NOTE: ALQ may drop below zero before we are able to meet the target + const std::string msg = fmt::format( + "Reduced (oil_rate, alq) from ({}, {}) to ({}, {}) to meet target " + "at {}. ", orig_oil_rate, orig_alq, oil_rate, alq, target); + displayDebugMessage_(msg); + } + else if (alq == orig_alq) { + // We might not be able to reduce ALQ, for example if ALQ starts out at zero. + const std::string msg = fmt::format("Not able to reduce ALQ {} further. " + "Oil rate is {} and oil target is {}", orig_alq, oil_rate, target); + 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 GasLiftWellState or nullptr +// +std::unique_ptr +GasLiftSingleWellGeneric:: +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 = false; + bool oil_is_limited = false; + bool gas_is_limited = false; + 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_; + double new_oil_rate, new_gas_rate, new_alq; + bool new_oil_is_limited, new_gas_is_limited; + std::tie(new_oil_rate, new_gas_rate, new_alq, new_oil_is_limited, new_gas_is_limited) + = maybeAdjustALQbeforeOptimizeLoop_( + increase, cur_alq, oil_rate, gas_rate, + oil_is_limited, gas_is_limited, potentials); + double delta_oil = 0.0; + double delta_gas = 0.0; + double delta_alq = 0.0; + OptimizeState state {*this, increase}; + if (checkInitialALQmodified_(new_alq, cur_alq)) { + delta_oil = new_oil_rate - oil_rate; + delta_gas = new_gas_rate - gas_rate; + delta_alq = new_alq - cur_alq; + if (!(state.checkGroupTargetsViolated(delta_oil, delta_gas)) && + !(state.checkGroupALQrateExceeded(delta_alq))) + { + oil_rate = new_oil_rate; + gas_rate = new_gas_rate; + oil_is_limited = new_oil_is_limited; + gas_is_limited = new_gas_is_limited; + cur_alq = new_alq; + success = true; + } + else { + state.stop_iteration = true; + } + } + auto temp_alq = cur_alq; + 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.checkGroupTargetsViolated(delta_oil, delta_gas)) 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; + delta_alq = *alq_opt - temp_alq; + if (state.checkGroupALQrateExceeded(delta_alq)) 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); + std::tie(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; + } + } +*/ + std::tie(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; + delta_oil = new_oil_rate - oil_rate; + delta_gas = new_gas_rate - gas_rate; + oil_rate = new_oil_rate; + gas_rate = new_gas_rate; + oil_is_limited = new_oil_is_limited; + gas_is_limited = new_gas_is_limited; + state.updateGroupRates(delta_oil, delta_gas, delta_alq); + } + 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; +} + +std::unique_ptr +GasLiftSingleWellGeneric:: +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; +} + +std::unique_ptr +GasLiftSingleWellGeneric:: +runOptimize2_() +{ + std::unique_ptr state; + state = tryIncreaseLiftGas_(); + if (!state || !(state->alqChanged())) { + state = tryDecreaseLiftGas_(); + } + return state; +} + +std::unique_ptr +GasLiftSingleWellGeneric:: +tryDecreaseLiftGas_() +{ + return runOptimizeLoop_(/*increase=*/ false); +} + +std::unique_ptr +GasLiftSingleWellGeneric:: +tryIncreaseLiftGas_() +{ + return runOptimizeLoop_(/*increase=*/ true); +} + +void +GasLiftSingleWellGeneric:: +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."); + } + } +} + +// Called when we should use a fixed ALQ value +void +GasLiftSingleWellGeneric:: +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. +bool +GasLiftSingleWellGeneric:: +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; + } +} void GasLiftSingleWellGeneric:: warnMaxIterationsExceeded_() @@ -702,6 +1092,70 @@ checkAlqOutsideLimits(double alq, [[maybe_unused]] double oil_rate) return result; } +bool +GasLiftSingleWellGeneric::OptimizeState:: +checkGroupALQrateExceeded(double delta_alq) +{ + const auto &pairs = + this->parent.group_info_.getWellGroups(this->parent.well_name_); + for (const auto &[group_name, efficiency] : pairs) { + auto max_alq_opt = this->parent.group_info_.maxAlq(group_name); + if (max_alq_opt) { + double alq = + this->parent.group_info_.alqRate(group_name) + efficiency * delta_alq; + if (alq > *max_alq_opt) { + if (this->parent.debug_) { + const std::string msg = fmt::format( + "Group {} : alq {} exceeds max_alq {}. Stopping iteration", + group_name, alq, *max_alq_opt); + this->parent.displayDebugMessage_(msg); + } + return true; + } + } + } + return false; +} + +bool +GasLiftSingleWellGeneric::OptimizeState:: +checkGroupTargetsViolated(double delta_oil, double delta_gas) +{ + const auto &pairs = + this->parent.group_info_.getWellGroups(this->parent.well_name_); + for (const auto &[group_name, efficiency] : pairs) { + auto oil_target_opt = this->parent.group_info_.oilTarget(group_name); + if (oil_target_opt) { + double oil_rate = + this->parent.group_info_.oilRate(group_name) + efficiency * delta_oil; + if (oil_rate > *oil_target_opt) { + if (this->parent.debug_) { + const std::string msg = fmt::format( + "Group {} : oil rate {} exceeds oil target {}. Stopping iteration", + group_name, oil_rate, *oil_target_opt); + this->parent.displayDebugMessage_(msg); + } + return true; + } + } + auto gas_target_opt = this->parent.group_info_.gasTarget(group_name); + if (gas_target_opt) { + double gas_rate = + this->parent.group_info_.gasRate(group_name) + efficiency * delta_gas; + if (gas_rate > *gas_target_opt) { + if (this->parent.debug_) { + const std::string msg = fmt::format( + "Group {} : gas rate {} exceeds gas target {}. Stopping iteration", + group_name, gas_rate, *gas_target_opt); + this->parent.displayDebugMessage_(msg); + } + return true; + } + } + } + return false; +} + bool GasLiftSingleWellGeneric::OptimizeState:: checkNegativeOilRate(double oil_rate) @@ -840,340 +1294,19 @@ getBhpWithLimit() return new_bhp; } -std::tuple -GasLiftSingleWellGeneric:: -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_) { - assert( alq <= orig_alq ); - if (alq < orig_alq) { - // NOTE: ALQ may drop below zero before we are able to meet the target - const std::string msg = fmt::format( - "Reduced (oil_rate, alq) from ({}, {}) to ({}, {}) to meet target " - "at {}. ", orig_oil_rate, orig_alq, oil_rate, alq, target); - displayDebugMessage_(msg); - } - else if (alq == orig_alq) { - // We might not be able to reduce ALQ, for example if ALQ starts out at zero. - const std::string msg = fmt::format("Not able to reduce ALQ {} further. " - "Oil rate is {} and oil target is {}", orig_alq, oil_rate, target); - 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 -// -std::unique_ptr -GasLiftSingleWellGeneric:: -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 = false; - bool oil_is_limited = false; - bool gas_is_limited = false; - 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; -} - -std::unique_ptr -GasLiftSingleWellGeneric:: -runOptimize(const int iteration_idx) -{ - 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, iteration_idx); - this->well_state_.setALQ(this->well_name_, alq); - } - } - } - return state; -} - -std::unique_ptr -GasLiftSingleWellGeneric:: -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; -} - -std::unique_ptr -GasLiftSingleWellGeneric:: -runOptimize2_() -{ - std::unique_ptr state; - state = tryIncreaseLiftGas_(); - if (!state || !(state->alqChanged())) { - state = tryDecreaseLiftGas_(); - } - return state; -} - -std::unique_ptr -GasLiftSingleWellGeneric:: -tryDecreaseLiftGas_() -{ - return runOptimizeLoop_(/*increase=*/ false); -} - -std::unique_ptr -GasLiftSingleWellGeneric:: -tryIncreaseLiftGas_() -{ - return runOptimizeLoop_(/*increase=*/ true); -} - void -GasLiftSingleWellGeneric:: -setAlqMinRate_(const GasLiftOpt::Well& well) +GasLiftSingleWellGeneric::OptimizeState:: +updateGroupRates(double delta_oil, double delta_gas, double delta_alq) { - // 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."); - } + const auto &pairs = + this->parent.group_info_.getWellGroups(this->parent.well_name_); + for (const auto &[group_name, efficiency] : pairs) { + int idx = this->parent.group_info_.getGroupIdx(group_name); + this->parent.sync_groups_.insert(idx); + this->parent.group_info_.update(group_name, + efficiency * delta_oil, efficiency * delta_gas, efficiency * delta_alq); } } -// Called when we should use a fixed ALQ value -void -GasLiftSingleWellGeneric:: -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. -bool -GasLiftSingleWellGeneric:: -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; - } -} } // namespace Opm diff --git a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp index 3b492e180..1bc6637b8 100644 --- a/opm/simulators/wells/GasLiftSingleWellGeneric.hpp +++ b/opm/simulators/wells/GasLiftSingleWellGeneric.hpp @@ -20,10 +20,14 @@ #ifndef OPM_GASLIFT_SINGLE_WELL_GENERIC_HEADER_INCLUDED #define OPM_GASLIFT_SINGLE_WELL_GENERIC_HEADER_INCLUDED +#include +#include + #include #include #include +#include #include #include @@ -51,6 +55,7 @@ protected: static constexpr double ALQ_EPSILON = 1e-8; public: + using GLiftSyncGroups = std::set; struct GradInfo { GradInfo() { } @@ -86,12 +91,16 @@ public: virtual const WellInterfaceGeneric& getStdWell() const = 0; protected: - GasLiftSingleWellGeneric(DeferredLogger &deferred_logger, - WellState &well_state, - const Well& ecl_well, - const SummaryState& summary_state, - const Schedule& schedule, - const int report_step_idx); + GasLiftSingleWellGeneric( + DeferredLogger &deferred_logger, + WellState &well_state, + const Well& ecl_well, + const SummaryState& summary_state, + GasLiftGroupInfo &group_info, + const Schedule& schedule, + const int report_step_idx, + GLiftSyncGroups &sync_groups + ); struct OptimizeState { @@ -114,6 +123,8 @@ protected: double gas_rate, double new_gas_rate); bool checkAlqOutsideLimits(double alq, double oil_rate); bool checkEcoGradient(double gradient); + bool checkGroupALQrateExceeded(double delta_alq); + bool checkGroupTargetsViolated(double delta_oil, double delta_gas); bool checkNegativeOilRate(double oil_rate); bool checkOilRateExceedsTarget(double oil_rate); bool checkRate(double rate, double limit, const std::string &rate_str) const; @@ -121,6 +132,7 @@ protected: bool computeBhpAtThpLimit(double alq); void debugShowIterationInfo(double alq); double getBhpWithLimit(); + void updateGroupRates(double delta_oil, double delta_gas, double delta_alq); void warn_(std::string msg) {parent.displayWarning_(msg);} }; @@ -183,15 +195,13 @@ protected: void logSuccess_(double alq, const int iteration_idx); - + std::tuple + maybeAdjustALQbeforeOptimizeLoop_( + bool increase, double alq, double oil_rate, double gas_rate, + bool oil_is_limited, bool gas_is_limited, std::vector &potentials); std::tuple - reduceALQtoOilTarget_(double alq, - double oil_rate, - double gas_rate, - bool oil_is_limited, - bool gas_is_limited, - std::vector& potentials); - + 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_(); @@ -212,7 +222,8 @@ protected: WellState& well_state_; const Well& ecl_well_; const SummaryState& summary_state_; - + GasLiftGroupInfo& group_info_; + GLiftSyncGroups& sync_groups_; const Well::ProductionControls controls_; double increment_; diff --git a/opm/simulators/wells/GasLiftSingleWell_impl.hpp b/opm/simulators/wells/GasLiftSingleWell_impl.hpp index bf0b75b3c..62209cca2 100644 --- a/opm/simulators/wells/GasLiftSingleWell_impl.hpp +++ b/opm/simulators/wells/GasLiftSingleWell_impl.hpp @@ -18,16 +18,6 @@ */ namespace Opm { -#include -#include -#include -#include -#include -#include -#include - -#include -#include template GasLiftSingleWell:: @@ -35,13 +25,22 @@ GasLiftSingleWell(const StdWell &std_well, const Simulator &ebos_simulator, const SummaryState &summary_state, DeferredLogger &deferred_logger, - WellState &well_state) - : GasLiftSingleWellGeneric(deferred_logger, - well_state, - std_well.wellEcl(), - summary_state, - ebos_simulator.vanguard().schedule(), - ebos_simulator.episodeIndex()) + WellState &well_state, + GasLiftGroupInfo &group_info, + GLiftSyncGroups &sync_groups + ) + // The parent class GasLiftSingleWellGeneric contains all stuff + // that is not dependent on TypeTag + : GasLiftSingleWellGeneric( + deferred_logger, + well_state, + std_well.wellEcl(), + summary_state, + group_info, + ebos_simulator.vanguard().schedule(), + ebos_simulator.episodeIndex(), + sync_groups + ) , ebos_simulator_{ebos_simulator} , std_well_{std_well} { diff --git a/opm/simulators/wells/GasLiftStage2.cpp b/opm/simulators/wells/GasLiftStage2.cpp index 1b357d6ab..6c7a147a4 100644 --- a/opm/simulators/wells/GasLiftStage2.cpp +++ b/opm/simulators/wells/GasLiftStage2.cpp @@ -458,7 +458,7 @@ mpiSyncGlobalGradVector_(std::vector &grads_global) const std::vector grads_local; for (auto itr = grads_global.begin(); itr != grads_global.end(); itr++) { - if (well_state_map_.count(itr->first) > 0) { + if (this->well_state_map_.count(itr->first) > 0) { grads_local.push_back(*itr); } } @@ -640,15 +640,18 @@ redistributeALQ_(std::vector &wells, const Group &group, std::vector &inc_grads, std::vector &dec_grads) { OptimizeState state {*this, group}; - // NOTE: 'inc_grads' and 'dec_grads' can never grow larger than wells.size() - // By reserving space here, we can ensure that any push_back() on these - // will never reallocate memory and invalidate any iterators. - inc_grads.reserve(wells.size()); - dec_grads.reserve(wells.size()); if (this->comm_.size() == 1) { + // NOTE: 'inc_grads' and 'dec_grads' can never grow larger than wells.size() + // By reserving space here, we can ensure that any push_back() on these + // will never reallocate memory and invalidate any iterators. + inc_grads.reserve(wells.size()); + dec_grads.reserve(wells.size()); state.calculateEcoGradients(wells, inc_grads, dec_grads); } else { + auto max_size = this->comm_.sum(wells.size()); + inc_grads.reserve(max_size); + dec_grads.reserve(max_size); std::vector inc_grads_local; std::vector dec_grads_local; inc_grads_local.reserve(wells.size()); @@ -662,7 +665,7 @@ redistributeALQ_(std::vector &wells, const Group &group, 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 + // need to calculate it (see above) since inc_grads and dec_grads are returned // and will be used by removeSurplusALQ_() later. return; } diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp index 7a96058c7..93b49c132 100644 --- a/opm/simulators/wells/GasLiftStage2.hpp +++ b/opm/simulators/wells/GasLiftStage2.hpp @@ -45,168 +45,172 @@ class Schedule; class WellInterfaceGeneric; class WellState; - class GasLiftStage2 { - using GasLiftSingleWell = GasLiftSingleWellGeneric; - using GLiftOptWells = std::map>; - using GLiftProdWells = std::map; - using GLiftWellStateMap = std::map>; - using GradPair = std::pair; - using GradPairItr = std::vector::iterator; - using GradInfo = typename GasLiftSingleWellGeneric::GradInfo; - using GradMap = std::map; - using MPIComm = typename Dune::MPIHelper::MPICommunicator; +class GasLiftStage2 { + using GasLiftSingleWell = GasLiftSingleWellGeneric; + using GLiftOptWells = std::map>; + using GLiftProdWells = std::map; + using GLiftWellStateMap = std::map>; + using GradPair = std::pair; + using GradPairItr = std::vector::iterator; + using GradInfo = typename GasLiftSingleWellGeneric::GradInfo; + using GradMap = std::map; + using MPIComm = typename Dune::MPIHelper::MPICommunicator; #if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7) - using Communication = Dune::Communication; + using Communication = Dune::Communication; #else - using Communication = Dune::CollectiveCommunication; + using Communication = Dune::CollectiveCommunication; #endif - static constexpr int Water = BlackoilPhases::Aqua; - static constexpr int Oil = BlackoilPhases::Liquid; - static constexpr int Gas = BlackoilPhases::Vapour; - public: - GasLiftStage2( - const int report_step_idx, - const Communication& comm, - const Schedule& schedule, - const SummaryState& summary_state, - 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_(GasLiftWellState &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 Group &group); - std::array getCurrentGroupRatesRecursive_( - const Group &group); - std::tuple getCurrentWellRates_( - const std::string &well_name, const std::string &group_name); - std::vector getGroupGliftWells_( - const Group &group); - void getGroupGliftWellsRecursive_( - const Group &group, std::vector &wells); - std::pair getStdWellRates_(const WellInterfaceGeneric &well); - void optimizeGroup_(const Group &group); - void optimizeGroupsRecursive_(const Group &group); - void recalculateGradientAndUpdateData_( - GradPairItr &grad_itr, bool increase, - std::vector &grads, std::vector &other_grads); - void redistributeALQ_( - std::vector &wells, const Group &group, - std::vector &inc_grads, std::vector &dec_grads); - void removeSurplusALQ_( - const Group &group, - std::vector &inc_grads, std::vector &dec_grads); - void saveGrad_(GradMap &map, const std::string &name, GradInfo &grad); - 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); - void mpiSyncGlobalGradVector_(std::vector &grads_global) const; - void mpiSyncLocalToGlobalGradVector_( - const std::vector &grads_local, - std::vector &grads_global) const; + static const int Water = BlackoilPhases::Aqua; + static const int Oil = BlackoilPhases::Liquid; + static const int Gas = BlackoilPhases::Vapour; +public: + GasLiftStage2( + const int report_step_idx, + const Communication& comm, + const Schedule& schedule, + const SummaryState& summary_state, + 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_(GasLiftWellState& 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 Group& group); + std::array getCurrentGroupRatesRecursive_( + const Group& group); + std::tuple getCurrentWellRates_( + const std::string& well_name, const std::string& group_name); + std::vector getGroupGliftWells_( + const Group& group); + void getGroupGliftWellsRecursive_( + const Group& group, std::vector& wells); + std::pair getStdWellRates_(const WellInterfaceGeneric& well); + void optimizeGroup_(const Group& group); + void optimizeGroupsRecursive_(const Group& group); + void recalculateGradientAndUpdateData_( + GradPairItr& grad_itr, bool increase, + std::vector& grads, std::vector& other_grads); + void redistributeALQ_( + std::vector& wells, const Group& group, + std::vector& inc_grads, std::vector& dec_grads); + void removeSurplusALQ_( + const Group& group, + std::vector& inc_grads, std::vector& dec_grads); + void saveGrad_(GradMap& map, const std::string& name, GradInfo& grad); + 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); + void mpiSyncGlobalGradVector_(std::vector& grads_global) const; + void mpiSyncLocalToGlobalGradVector_( + const std::vector& grads_local, + std::vector& grads_global) const; - DeferredLogger &deferred_logger_; - WellState &well_state_; - GLiftProdWells &prod_wells_; - GLiftOptWells &stage1_wells_; - GLiftWellStateMap &well_state_map_; + DeferredLogger& deferred_logger_; + 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 GasLiftOpt& glo_; - const Communication &comm_; - GradMap inc_grads_; - GradMap dec_grads_; - bool debug_; - static constexpr int max_iterations_ = 1000; - //int time_step_idx_; + int report_step_idx_; + const SummaryState& summary_state_; + const Schedule& schedule_; + const GasLiftOpt& glo_; + const Communication& comm_; + GradMap inc_grads_; + GradMap dec_grads_; + bool debug_; + int max_iterations_ = 1000; + //int time_step_idx_; - struct OptimizeState { - OptimizeState( GasLiftStage2 &parent_, const Group &group_ ) : - parent{parent_}, - group{group_}, - it{0} - {} - GasLiftStage2 &parent; - const Group &group; - int it; + struct OptimizeState { + OptimizeState(GasLiftStage2& parent_, const Group& group_ ) : + parent{parent_}, + group{group_}, + it{0} + {} + GasLiftStage2& parent; + const Group& group; + int it; - 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); + 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); + private: + void displayDebugMessage_(const std::string& msg); + void displayWarning_(const std::string& msg); - }; + }; - struct SurplusState { - SurplusState( GasLiftStage2 &parent_, const 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 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; + struct SurplusState { + SurplusState( GasLiftStage2& parent_, const 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 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); - }; - }; + 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); + }; +}; } // namespace Opm diff --git a/opm/simulators/wells/MultisegmentWell.hpp b/opm/simulators/wells/MultisegmentWell.hpp index 17ad22de3..89f976aab 100644 --- a/opm/simulators/wells/MultisegmentWell.hpp +++ b/opm/simulators/wells/MultisegmentWell.hpp @@ -56,6 +56,7 @@ namespace Opm using typename Base::GLiftProdWells; using typename Base::GLiftOptWells; using typename Base::GLiftWellStateMap; + using typename Base::GLiftSyncGroups; /// the number of reservior equations using Base::numEq; @@ -108,7 +109,9 @@ namespace Opm DeferredLogger&, GLiftProdWells &, GLiftOptWells &, - GLiftWellStateMap & + GLiftWellStateMap &, + GasLiftGroupInfo &, + GLiftSyncGroups & ) const override { // Not implemented yet } diff --git a/opm/simulators/wells/StandardWell.hpp b/opm/simulators/wells/StandardWell.hpp index aaf0d6321..012bd6c64 100644 --- a/opm/simulators/wells/StandardWell.hpp +++ b/opm/simulators/wells/StandardWell.hpp @@ -32,6 +32,7 @@ #include #include #include +#include #include #include @@ -85,6 +86,7 @@ namespace Opm using typename Base::GLiftOptWells; using typename Base::GLiftProdWells; using typename Base::GLiftWellStateMap; + using typename Base::GLiftSyncGroups; using Base::numEq; using Base::numPhases; @@ -230,7 +232,9 @@ namespace Opm DeferredLogger& deferred_logger, GLiftProdWells &prod_wells, GLiftOptWells &glift_wells, - GLiftWellStateMap &state_map + GLiftWellStateMap &state_map, + GasLiftGroupInfo &group_info, + GLiftSyncGroups &sync_groups ) const override; /* returns BHP */ diff --git a/opm/simulators/wells/StandardWell_impl.hpp b/opm/simulators/wells/StandardWell_impl.hpp index 60ab48a59..c3742eda1 100644 --- a/opm/simulators/wells/StandardWell_impl.hpp +++ b/opm/simulators/wells/StandardWell_impl.hpp @@ -1499,34 +1499,24 @@ namespace Opm DeferredLogger& deferred_logger, GLiftProdWells &prod_wells, GLiftOptWells &glift_wells, - GLiftWellStateMap &glift_state_map - //std::map &prod_wells + GLiftWellStateMap &glift_state_map, + GasLiftGroupInfo &group_info, + GLiftSyncGroups &sync_groups ) const { - const auto& well = well_ecl_; - if (well.isProducer()) { - const auto& summary_state = ebos_simulator.vanguard().summaryState(); - if ( this->Base::wellHasTHPConstraints(summary_state) ) { - const int iteration_idx = ebos_simulator.model().newtonMethod().numIterations(); - if (this->doGasLiftOptimize(well_state, - ebos_simulator.episodeIndex(), - iteration_idx, - ebos_simulator.vanguard().schedule(), - deferred_logger)) { - std::unique_ptr glift - = std::make_unique( - *this, ebos_simulator, summary_state, - deferred_logger, well_state); - auto state = glift->runOptimize(iteration_idx); - 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}); + const auto& summary_state = ebos_simulator.vanguard().summaryState(); + std::unique_ptr glift + = std::make_unique( + *this, ebos_simulator, summary_state, + deferred_logger, well_state, group_info, sync_groups); + auto state = glift->runOptimize( + ebos_simulator.model().newtonMethod().numIterations()); + 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}); } template diff --git a/opm/simulators/wells/WellInterface.hpp b/opm/simulators/wells/WellInterface.hpp index 05229a599..9d7e11b6d 100644 --- a/opm/simulators/wells/WellInterface.hpp +++ b/opm/simulators/wells/WellInterface.hpp @@ -42,7 +42,9 @@ namespace Opm { template class GasLiftSingleWell; template class BlackoilWellModel; } +#include #include +#include #include #include @@ -84,6 +86,7 @@ public: using GLiftProdWells = typename BlackoilWellModel::GLiftProdWells; using GLiftWellStateMap = typename BlackoilWellModel::GLiftWellStateMap; + using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; static const int numEq = Indices::numEq; static const int numPhases = Indices::numPhases; @@ -167,7 +170,9 @@ public: DeferredLogger& deferred_logger, GLiftProdWells& prod_wells, GLiftOptWells& glift_wells, - GLiftWellStateMap& state_map + GLiftWellStateMap& state_map, + GasLiftGroupInfo &group_info, + GLiftSyncGroups &sync_groups ) const = 0; /// using the solution x to recover the solution xw for wells and applying diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index c52216b6e..6a44ad9be 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -250,6 +250,8 @@ void WellState::init(const std::vector& cellPressures, const int nw = wells_ecl.size(); + do_glift_optimization_ = true; + if( nw == 0 ) return ; // Initialize perfphaserates_, which must be done here. @@ -437,7 +439,6 @@ void WellState::init(const std::vector& cellPressures, updateWellsDefaultALQ(wells_ecl); - do_glift_optimization_ = true; } void WellState::resize(const std::vector& wells_ecl, diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 2acc2cf31..c476fd4c8 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -256,15 +256,6 @@ public: void gliftTimeStepInit() { this->alq_state.reset_count(); - disableGliftOptimization(); - } - - void disableGliftOptimization() { - do_glift_optimization_ = false; - } - - void enableGliftOptimization() { - do_glift_optimization_ = true; } int wellNameToGlobalIdx(const std::string &name) { diff --git a/tests/test_glift1.cpp b/tests/test_glift1.cpp index d4235b472..1e499f762 100644 --- a/tests/test_glift1.cpp +++ b/tests/test_glift1.cpp @@ -38,6 +38,8 @@ #include #include #include +#include +#include //#include //#include #include @@ -122,7 +124,11 @@ BOOST_AUTO_TEST_CASE(G1) using WellState = Opm::WellState; using StdWell = Opm::StandardWell; using GasLiftSingleWell = Opm::GasLiftSingleWell; + using GasLiftGroupInfo = Opm::GasLiftGroupInfo; + using GasLiftSingleWellGeneric = Opm::GasLiftSingleWellGeneric; + using GLiftEclWells = typename GasLiftGroupInfo::GLiftEclWells; const std::string filename = "GLIFT1.DATA"; + using GLiftSyncGroups = typename GasLiftSingleWellGeneric::GLiftSyncGroups; auto simulator = initSimulator(filename.data()); @@ -161,9 +167,23 @@ BOOST_AUTO_TEST_CASE(G1) BOOST_CHECK_EQUAL( well.name(), "B-1H"); const auto& summary_state = simulator->vanguard().summaryState(); WellState &well_state = well_model.wellState(); + GLiftEclWells ecl_well_map; + well_model.initGliftEclWellMap(ecl_well_map); + const int iteration_idx = simulator->model().newtonMethod().numIterations(); + GasLiftGroupInfo group_info { + ecl_well_map, + schedule, + summary_state, + simulator->episodeIndex(), + iteration_idx, + well_model.phaseUsage(), + deferred_logger, + well_state + }; + GLiftSyncGroups sync_groups; GasLiftSingleWell glift {*std_well, *(simulator.get()), summary_state, - deferred_logger, well_state}; - auto state = glift.runOptimize(simulator->model().newtonMethod().numIterations()); + deferred_logger, well_state, group_info, sync_groups}; + auto state = glift.runOptimize(iteration_idx); BOOST_CHECK_CLOSE(state->oilRate(), 0.01736111111111111, 1e-8); BOOST_CHECK_CLOSE(state->gasRate(), 1.6464646999768586, 1e-8); BOOST_CHECK(state->oilIsLimited()); @@ -172,3 +192,4 @@ BOOST_AUTO_TEST_CASE(G1) BOOST_CHECK_CLOSE(state->alq(), 0.0, 1e-8); BOOST_CHECK(!state->increase().has_value()); } +