diff --git a/opm/simulators/wells/GasLiftStage2.hpp b/opm/simulators/wells/GasLiftStage2.hpp index 9151f336b..fceb24529 100644 --- a/opm/simulators/wells/GasLiftStage2.hpp +++ b/opm/simulators/wells/GasLiftStage2.hpp @@ -69,6 +69,12 @@ namespace Opm using GradPairItr = std::vector::iterator; using GradInfo = typename GasLiftSingleWell::GradInfo; using GradMap = std::map; + using MPIComm = typename Dune::MPIHelper::MPICommunicator; +#if DUNE_VERSION_NEWER(DUNE_COMMON, 2, 7) + using Communication = Dune::Communication; +#else + using Communication = Dune::CollectiveCommunication; +#endif static const int Water = BlackoilPhases::Aqua; static const int Oil = BlackoilPhases::Liquid; static const int Gas = BlackoilPhases::Vapour; @@ -127,8 +133,10 @@ namespace Opm const std::string &name, GradInfo &grad, bool increase); void updateGradVector_( const std::string &name, std::vector &grads, double grad); - std::vector updateGlobalGradVector_(const std::vector &grads_global) const; - std::vector localToGlobalGradVector_(const std::vector &grads_local) const; + void mpiSyncGlobalGradVector_(std::vector &grads_global) const; + void mpiSyncLocalToGlobalGradVector_( + const std::vector &grads_local, + std::vector &grads_global) const; DeferredLogger &deferred_logger_; @@ -144,6 +152,7 @@ namespace Opm const Schedule &schedule_; const PhaseUsage &phase_usage_; const GasLiftOpt& glo_; + const Communication &comm_; GradMap inc_grads_; GradMap dec_grads_; bool debug_; diff --git a/opm/simulators/wells/GasLiftStage2_impl.hpp b/opm/simulators/wells/GasLiftStage2_impl.hpp index 9b8d75c73..630c35e20 100644 --- a/opm/simulators/wells/GasLiftStage2_impl.hpp +++ b/opm/simulators/wells/GasLiftStage2_impl.hpp @@ -56,6 +56,7 @@ GasLiftStage2( schedule_{ebos_simulator.vanguard().schedule()}, phase_usage_{well_model_.phaseUsage()}, glo_{schedule_.glo(report_step_idx_)}, + comm_{ebos_simulator.vanguard().grid().comm()}, debug_{false} { // this->time_step_idx_ @@ -286,8 +287,7 @@ GasLiftStage2:: getCurrentGroupRates_(const Group &group) { auto rates = getCurrentGroupRatesRecursive_(group); - const auto& comm = ebos_simulator_.vanguard().grid().comm(); - comm.sum(rates.data(), rates.size()); + this->comm_.sum(rates.data(), rates.size()); auto [oil_rate, gas_rate, alq] = rates; if (this->debug_) { const std::string msg = fmt::format( @@ -468,6 +468,64 @@ getGroupGliftWellsRecursive_(const Group &group, } } +template +void +GasLiftStage2:: +mpiSyncGlobalGradVector_(std::vector &grads_global) const +{ + if (this->comm_.size() == 1) + return; + + std::vector grads_local; + for (auto itr = grads_global.begin(); itr != grads_global.end(); itr++) { + if (well_state_map_.count(itr->first) > 0) { + grads_local.push_back(*itr); + } + } + mpiSyncLocalToGlobalGradVector_(grads_local, grads_global); +} + +template +void +GasLiftStage2:: +mpiSyncLocalToGlobalGradVector_( + const std::vector &grads_local, std::vector &grads_global) const +{ + assert(this->comm_.size() > 1); // The parent should check if comm. size is > 1 + using Pair = std::pair; + std::vector grads_local_tmp; + grads_local_tmp.reserve(grads_local.size()); + for (size_t i = 0; i < grads_local.size(); ++i) { + if(!this->well_state_.wellIsOwned(grads_local[i].first)) + continue; + grads_local_tmp.push_back( + std::make_pair( + this->well_state_.wellNameToGlobalIdx(grads_local[i].first), + grads_local[i].second)); + } + + std::vector sizes_(this->comm_.size()); + std::vector displ_(this->comm_.size() + 1, 0); + int mySize = grads_local_tmp.size(); + this->comm_.allgather(&mySize, 1, sizes_.data()); + std::partial_sum(sizes_.begin(), sizes_.end(), displ_.begin()+1); + std::vector grads_global_tmp(displ_.back()); + + this->comm_.allgatherv(grads_local_tmp.data(), grads_local_tmp.size(), + grads_global_tmp.data(), sizes_.data(), displ_.data()); + + // NOTE: This leaves the capacity of 'grads_global' unchanged, so + // memory is not reallocated here + grads_global.clear(); + + for (size_t i = 0; i < grads_global_tmp.size(); ++i) { + grads_global.emplace_back( + std::make_pair( + well_state_.globalIdxToWellName(grads_global_tmp[i].first), + grads_global_tmp[i].second)); + } +} + template void GasLiftStage2:: @@ -514,13 +572,16 @@ optimizeGroupsRecursive_(const Group &group) optimizeGroup_(group); } - template void GasLiftStage2:: recalculateGradientAndUpdateData_( GradPairItr &grad_itr, bool increase, - std::vector &grads, std::vector &other_grads) + + //incremental and decremental gradients, if 'grads' are incremental, then + // 'other_grads' are decremental, or conversely, if 'grads' are decremental, then + // 'other_grads' are incremental + std::vector &grads, std::vector &other_grads) { // NOTE: We make a copy of the name string instead of taking a reference // since we may have to erase grad_itr (in the "else" condition below) @@ -547,6 +608,8 @@ recalculateGradientAndUpdateData_( // The old incremental gradient becomes the new decremental gradient // or the old decremental gradient becomes the new incremental gradient updateGrad_(name, *old_grad, !increase); + // NOTE: This may invalidate any iterator into 'other_grads' since + // updateGradVector_() will do a push_back() if 'name' is not found.. updateGradVector_(name, other_grads, old_grad->grad); } } @@ -601,15 +664,24 @@ redistributeALQ_(std::vector &wells, const Group &group, std::vector &inc_grads, std::vector &dec_grads) { OptimizeState state {*this, group}; - std::vector inc_grads_local; - std::vector dec_grads_local; - inc_grads_local.reserve(wells.size()); - dec_grads_local.reserve(wells.size()); - state.calculateEcoGradients(wells, inc_grads_local, dec_grads_local); - - // the gradients needs to be communicated to all ranks - dec_grads = localToGlobalGradVector_(dec_grads_local); - inc_grads = localToGlobalGradVector_(inc_grads_local); + // 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) { + state.calculateEcoGradients(wells, inc_grads, dec_grads); + } + else { + std::vector inc_grads_local; + std::vector dec_grads_local; + inc_grads_local.reserve(wells.size()); + dec_grads_local.reserve(wells.size()); + state.calculateEcoGradients(wells, inc_grads_local, dec_grads_local); + // the gradients needs to be communicated to all ranks + mpiSyncLocalToGlobalGradVector_(dec_grads_local, dec_grads); + mpiSyncLocalToGlobalGradVector_(inc_grads_local, inc_grads); + } if (!state.checkAtLeastTwoWells(wells)) { // NOTE: Even though we here in redistributeALQ_() do not use the @@ -705,8 +777,8 @@ removeSurplusALQ_(const Group &group, dec_grad_itr, /*increase=*/false, dec_grads, inc_grads); // The dec_grads and inc_grads needs to be syncronized across ranks - dec_grads = updateGlobalGradVector_(dec_grads); - inc_grads = updateGlobalGradVector_(inc_grads); + mpiSyncGlobalGradVector_(dec_grads); + mpiSyncGlobalGradVector_(inc_grads); // NOTE: recalculateGradientAndUpdateData_() will remove the current gradient // from dec_grads if it cannot calculate a new decremental gradient. // This will invalidate dec_grad_itr and well_name @@ -803,58 +875,6 @@ updateGradVector_(const std::string &name, std::vector &grads, double // later in getEcoGradients() } -template -std::vector::GradPair> -GasLiftStage2:: -localToGlobalGradVector_(const std::vector &grads_local) const -{ - const auto& comm = ebos_simulator_.vanguard().grid().comm(); - if (comm.size() == 1) - return grads_local; - - using Pair = std::pair; - std::vector grads_local_tmp; - grads_local_tmp.reserve(grads_local.size()); - for (size_t i = 0; i < grads_local.size(); ++i) { - if(!well_state_.wellIsOwned(grads_local[i].first)) - continue; - - grads_local_tmp.push_back(std::make_pair(well_state_.wellNameToGlobalIdx(grads_local[i].first), grads_local[i].second)); - } - - std::vector sizes_(comm.size()); - std::vector displ_(comm.size() + 1, 0); - int mySize = grads_local_tmp.size(); - comm.allgather(&mySize, 1, sizes_.data()); - std::partial_sum(sizes_.begin(), sizes_.end(), displ_.begin()+1); - std::vector grads_global(displ_.back()); - - comm.allgatherv(grads_local_tmp.data(), grads_local_tmp.size(), grads_global.data(), sizes_.data(), displ_.data()); - std::vector grads(grads_global.size()); - for (size_t i = 0; i < grads_global.size(); ++i) { - grads[i] = std::make_pair(well_state_.globalIdxToWellName(grads_global[i].first), grads_global[i].second); - } - return grads; -} - -template -std::vector::GradPair> -GasLiftStage2:: -updateGlobalGradVector_(const std::vector &grads_global) const -{ - const auto& comm = ebos_simulator_.vanguard().grid().comm(); - if (comm.size() == 1) - return grads_global; - - std::vector grads_local; - for (auto itr = grads_global.begin(); itr != grads_global.end(); itr++) { - if (well_state_map_.count(itr->first) > 0) { - grads_local.push_back(*itr); - } - } - return localToGlobalGradVector_(grads_local); -} - /*********************************************** * Public methods declared in OptimizeState ***********************************************/ @@ -894,8 +914,7 @@ checkAtLeastTwoWells(std::vector &wells) continue; numberOfwells++; } - const auto& comm = this->parent.ebos_simulator_.vanguard().grid().comm(); - numberOfwells = comm.sum(numberOfwells); + numberOfwells = this->parent.comm_.sum(numberOfwells); if (numberOfwells < 2) { const std::string msg = fmt::format( "skipping: too few wells ({}) to optimize.", numberOfwells); @@ -970,8 +989,8 @@ recalculateGradients( min_dec_grad_itr, /*increase=*/false, dec_grads, inc_grads); // The dec_grads and inc_grads needs to be syncronized across ranks - dec_grads = this->parent.updateGlobalGradVector_(dec_grads); - inc_grads = this->parent.updateGlobalGradVector_(inc_grads); + this->parent.mpiSyncGlobalGradVector_(dec_grads); + this->parent.mpiSyncGlobalGradVector_(inc_grads); } // Take one ALQ increment from well1, and give it to well2 @@ -1132,8 +1151,7 @@ updateRates(const std::string &well_name) } // and communicate the results - const auto& comm = this->parent.ebos_simulator_.vanguard().grid().comm(); - comm.sum(delta.data(), delta.size()); + this->parent.comm_.sum(delta.data(), delta.size()); // and update const auto& [delta_oil, delta_gas, delta_alq] = delta;