Merge pull request #3148 from totto82/gliftmpi2

make gaslift parallel

opm/simulators/wells/GasLiftStage2.hpp

@@ -99,7 +99,7 @@ namespace Opm
         void displayWarning_(const std::string &msg);
         std::tuple<double, double, double> getCurrentGroupRates_(
             const Opm::Group &group);
-        std::tuple<double, double, double> getCurrentGroupRatesRecursive_(
+        std::array<double,3> getCurrentGroupRatesRecursive_(
             const Opm::Group &group);
         std::tuple<double, double, double> getCurrentWellRates_(
             const std::string &well_name, const std::string &group_name);
@@ -127,6 +127,9 @@ namespace Opm
             const std::string &name, GradInfo &grad, bool increase);
         void updateGradVector_(
             const std::string &name, std::vector<GradPair> &grads, double grad);
+        std::vector<GradPair> updateGlobalGradVector_(const std::vector<GradPair> &grads_global) const;
+        std::vector<GradPair> localToGlobalGradVector_(const std::vector<GradPair> &grads_local) const;
+
 
         DeferredLogger &deferred_logger_;
         const Simulator &ebos_simulator_;
@@ -213,7 +216,7 @@ namespace Opm
             bool checkEcoGradient(const std::string &well_name, double eco_grad);
             bool checkGasTarget();
             bool checkOilTarget();
-            void updateRates(const std::string &name, const GradInfo &gi);
+            void updateRates(const std::string &name);
         private:
         };
     };

opm/simulators/wells/GasLiftStage2_impl.hpp

@@ -100,8 +100,13 @@ void
 GasLiftStage2<TypeTag>::
 addOrRemoveALQincrement_(GradMap &grad_map, const std::string well_name, bool add)
 {
+    // only applies to wells in the well_state_map (i.e. wells on this rank)
+    auto it = this->well_state_map_.find(well_name);
+    if (it == this->well_state_map_.end())
+        return;
+
+    GLiftWellState &state = *(it->second.get());
     const GradInfo &gi = grad_map.at(well_name);
-    GLiftWellState &state = *(this->well_state_map_.at(well_name).get());
     if (this->debug_) {
         auto new_alq = gi.alq;
         auto old_alq = state.alq();
@@ -122,6 +127,11 @@ GasLiftStage2<TypeTag>::
 calcIncOrDecGrad_(
     const std::string well_name, const GasLiftSingleWell &gs_well, bool increase)
 {
+
+    // only applies to wells in the well_state_map (i.e. wells on this rank)
+    if(this->well_state_map_.count(well_name) == 0)
+        return std::nullopt;
+
     GLiftWellState &state = *(this->well_state_map_.at(well_name).get());
     if (checkRateAlreadyLimited_(state, increase)) {
         /*
@@ -273,23 +283,23 @@ std::tuple<double, double, double>
 GasLiftStage2<TypeTag>::
 getCurrentGroupRates_(const Opm::Group &group)
 {
+    auto rates = getCurrentGroupRatesRecursive_(group);
+    const auto& comm = ebos_simulator_.vanguard().grid().comm();
+    comm.sum(rates.data(), rates.size());
+    auto [oil_rate, gas_rate, alq] = rates;
     if (this->debug_) {
-        auto [oil_rate, gas_rate, alq] = getCurrentGroupRatesRecursive_(group);
         const std::string msg = fmt::format(
             "Current group rates for {} : oil: {}, gas: {}, alq: {}",
             group.name(), oil_rate, gas_rate, alq);
         displayDebugMessage2B_(msg);
-        return {oil_rate, gas_rate, alq};
-    }
-    else {
-        return getCurrentGroupRatesRecursive_(group);
     }
+
+    return {oil_rate, gas_rate, alq};
 }
 
 
 template<typename TypeTag>
-std::tuple<double, double, double>
-GasLiftStage2<TypeTag>::
+std::array <double, 3> GasLiftStage2<TypeTag>::
 getCurrentGroupRatesRecursive_(const Opm::Group &group)
 {
     double oil_rate = 0.0;
@@ -304,6 +314,7 @@ getCurrentGroupRatesRecursive_(const Opm::Group &group)
             gas_rate += sw_gas_rate;
             alq += sw_alq;
         }
+
     }
     else {
         for (const std::string& group_name : group.groups()) {
@@ -317,8 +328,8 @@ getCurrentGroupRatesRecursive_(const Opm::Group &group)
                 // they are added to the lift gas supply rate of the
                 // parent group.
                 const auto gefac = sub_group.getGroupEfficiencyFactor();
-                auto [sg_oil_rate, sg_gas_rate, sg_alq] =
-                    getCurrentGroupRatesRecursive_(sub_group);
+                auto rates = getCurrentGroupRatesRecursive_(sub_group);
+                auto [sg_oil_rate, sg_gas_rate, sg_alq] = rates;
                 oil_rate += (gefac * sg_oil_rate);
                 gas_rate += (gefac * sg_gas_rate);
                 alq += (gefac * sg_alq);
@@ -352,6 +363,14 @@ getCurrentWellRates_(const std::string &well_name, const std::string &group_name
         success = true;
         if ( this->debug_) debug_info = "(B)";
     }
+
+    if (well_ptr) {
+        // we only want rates from wells owned by the rank
+        if (!well_state_.wellIsOwned(well_ptr->indexOfWell(), well_name)) {
+            success = false;
+        }
+    }
+
     if (success) {
         assert(well_ptr);
         assert(well_ptr->isProducer());
@@ -504,17 +523,23 @@ recalculateGradientAndUpdateData_(
     // NOTE: We make a copy of the name string instead of taking a reference
     //   since we may have to erase grad_itr (in the "else" condition below)
     const std::string name = grad_itr->first;
-    GasLiftSingleWell &gs_well = *(this->stage1_wells_.at(name).get());
-    auto grad = calcIncOrDecGrad_(name, gs_well, increase);
     std::optional<GradInfo> old_grad = std::nullopt;
-    if (grad) {
-        grad_itr->second = grad->grad;
-        old_grad = updateGrad_(name, *grad, increase);
-    }
-    else {
-        grads.erase(grad_itr); // NOTE: this invalidates grad_itr
-        old_grad = deleteGrad_(name, increase);
+
+    // only applies to wells in the well_state_map (i.e. wells on this rank)
+    // the grads and other grads are synchronized later
+    if(this->stage1_wells_.count(name) > 0) {
+        GasLiftSingleWell &gs_well = *(this->stage1_wells_.at(name).get());
+        auto grad = calcIncOrDecGrad_(name, gs_well, increase);
+        if (grad) {
+            grad_itr->second = grad->grad;
+            old_grad = updateGrad_(name, *grad, increase);
+        }
+        else {
+            grads.erase(grad_itr); // NOTE: this invalidates grad_itr
+            old_grad = deleteGrad_(name, increase);
+        }
     }
+
     if (old_grad) {
         // NOTE: Either creates a new item or reassigns
         // The old incremental gradient becomes the new decremental gradient
@@ -574,9 +599,16 @@ redistributeALQ_(std::vector<GasLiftSingleWell *> &wells,  const Opm::Group &gro
     std::vector<GradPair> &inc_grads, std::vector<GradPair> &dec_grads)
 {
     OptimizeState state {*this, group};
-    inc_grads.reserve(wells.size());
-    dec_grads.reserve(wells.size());
-    state.calculateEcoGradients(wells, inc_grads, dec_grads);
+    std::vector<GradPair> inc_grads_local;
+    std::vector<GradPair> 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);
+
     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
@@ -584,7 +616,6 @@ redistributeALQ_(std::vector<GasLiftSingleWell *> &wells,  const Opm::Group &gro
         //   and will be used by removeSurplusALQ_() later.
         return;
     }
-    assert(wells.size() >= 2);
     bool stop_iteration = false;
     while (!stop_iteration && (state.it++ <= this->max_iterations_)) {
         state.debugShowIterationInfo();
@@ -645,12 +676,13 @@ removeSurplusALQ_(const Opm::Group &group,
     }
     SurplusState state {*this, group, oil_rate, gas_rate, alq,
             min_eco_grad, controls.oil_target, controls.gas_target, max_glift };
+
     while (!stop_iteration) {
         if (dec_grads.size() >= 2) {
             sortGradients_(dec_grads);
         }
         auto dec_grad_itr = dec_grads.begin();
-        const auto &well_name = dec_grad_itr->first;
+        const auto well_name = dec_grad_itr->first;
         auto eco_grad = dec_grad_itr->second;
         bool remove = false;
         if (state.checkOilTarget() || state.checkGasTarget() || state.checkALQlimit()) {
@@ -665,11 +697,14 @@ removeSurplusALQ_(const Opm::Group &group,
             if (state.checkEcoGradient(well_name, eco_grad)) remove = true;
         }
         if (remove) {
-            const GradInfo &gi = this->dec_grads_.at(well_name);
-            state.updateRates(well_name, gi);
+            state.updateRates(well_name);
             state.addOrRemoveALQincrement( this->dec_grads_, well_name, /*add=*/false);
             recalculateGradientAndUpdateData_(
-                dec_grad_itr, /*increase=*/false, dec_grads, inc_grads);
+                        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);
             // 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
@@ -766,6 +801,57 @@ updateGradVector_(const std::string &name, std::vector<GradPair> &grads, double
     //   later in getEcoGradients()
 }
 
+template<typename TypeTag>
+std::vector<typename GasLiftStage2<TypeTag>::GradPair>
+GasLiftStage2<TypeTag>::
+localToGlobalGradVector_(const std::vector<GradPair> &grads_local) const
+{
+    const auto& comm = ebos_simulator_.vanguard().grid().comm();
+    if (comm.size() == 1)
+        return grads_local;
+
+    using Pair = std::pair<int, double>;
+    std::vector<Pair> 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<int> sizes_(comm.size());
+    std::vector<int> 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<Pair> grads_global(displ_.back());
+
+    comm.allgatherv(grads_local_tmp.data(), grads_local_tmp.size(), grads_global.data(), sizes_.data(), displ_.data());
+    std::vector<GradPair> 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<typename TypeTag>
+std::vector<typename GasLiftStage2<TypeTag>::GradPair>
+GasLiftStage2<TypeTag>::
+updateGlobalGradVector_(const std::vector<GradPair> &grads_global) const
+{
+    const auto& comm = ebos_simulator_.vanguard().grid().comm();
+    if (comm.size() == 1)
+        return grads_global;
+
+    std::vector<GradPair> 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
@@ -799,9 +885,18 @@ bool
 GasLiftStage2<TypeTag>::OptimizeState::
 checkAtLeastTwoWells(std::vector<GasLiftSingleWell *> &wells)
 {
-    if (wells.size() < 2) {
+    int numberOfwells = 0;
+    for (auto well : wells){
+        int index_of_wells = well->getStdWell().indexOfWell();
+        if (!this->parent.well_state_.wellIsOwned(index_of_wells, well->name()))
+            continue;
+        numberOfwells++;
+    }
+    const auto& comm = this->parent.ebos_simulator_.vanguard().grid().comm();
+    numberOfwells = comm.sum(numberOfwells);
+    if (numberOfwells < 2) {
         const std::string msg = fmt::format(
-            "skipping: too few wells ({}) to optimize.", wells.size());
+            "skipping: too few wells ({}) to optimize.", numberOfwells);
         displayDebugMessage_(msg);
         return false;
     }
@@ -871,6 +966,10 @@ recalculateGradients(
         max_inc_grad_itr, /*increase=*/true, inc_grads, dec_grads);
     this->parent.recalculateGradientAndUpdateData_(
         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);
 }
 
 // Take one ALQ increment from well1, and give it to well2
@@ -1010,17 +1109,35 @@ checkOilTarget()
 template<typename TypeTag>
 void
 GasLiftStage2<TypeTag>::SurplusState::
-updateRates(const std::string &well_name, const GradInfo &gi)
+updateRates(const std::string &well_name)
 {
-    GLiftWellState &state = *(this->parent.well_state_map_.at(well_name).get());
-    GasLiftSingleWell &gs_well = *(this->parent.stage1_wells_.at(well_name).get());
-    const WellInterface<TypeTag> &well = gs_well.getStdWell();
-    const auto &well_ecl = well.wellEcl();
-    double factor = well_ecl.getEfficiencyFactor();
-    auto delta_oil = gi.new_oil_rate - state.oilRate();
-    this->oil_rate += factor * delta_oil;
-    auto delta_gas = gi.new_gas_rate - state.gasRate();
-    this->gas_rate += factor * delta_gas;
-    auto delta_alq = gi.alq - state.alq();
-    this->alq += factor * delta_alq;
+    std::array<double, 3> delta = {0.0,0.0,0.0};
+    // compute the delta on wells on own rank
+    if (this->parent.well_state_map_.count(well_name) > 0) {
+        const GradInfo &gi = this->parent.dec_grads_.at(well_name);
+        GLiftWellState &state = *(this->parent.well_state_map_.at(well_name).get());
+        GasLiftSingleWell &gs_well = *(this->parent.stage1_wells_.at(well_name).get());
+        const WellInterface<TypeTag> &well = gs_well.getStdWell();
+        // only get deltas for wells owned by this rank
+        if (this->parent.well_state_.wellIsOwned(well.indexOfWell(), well_name)) {
+            const auto &well_ecl = well.wellEcl();
+            double factor = well_ecl.getEfficiencyFactor();
+            auto& [delta_oil, delta_gas, delta_alq] = delta;
+                    delta_oil = factor * (gi.new_oil_rate - state.oilRate());
+                    delta_gas = factor * (gi.new_gas_rate - state.gasRate());
+                    delta_alq = factor * (gi.alq - state.alq());
+        }
+    }
+
+    // and communicate the results
+    const auto& comm = this->parent.ebos_simulator_.vanguard().grid().comm();
+    comm.sum(delta.data(), delta.size());
+
+    // and update
+    const auto& [delta_oil, delta_gas, delta_alq] = delta;
+    this->oil_rate += delta_oil;
+    this->gas_rate += delta_gas;
+    this->alq += delta_alq;
 }
+
+

opm/simulators/wells/StandardWell_impl.hpp

@@ -2692,13 +2692,6 @@ namespace Opm
     {
 
         gliftDebug("checking if GLIFT should be done..", deferred_logger);
-        std::size_t num_procs = ebos_simulator.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 (!well_state.gliftOptimizationEnabled()) {
             gliftDebug("Optimization disabled in WellState", deferred_logger);
             return false;

opm/simulators/wells/WellState.hpp

@@ -209,6 +209,16 @@ namespace Opm
             return well_info.isOwner();
         }
 
+        bool wellIsOwned(const std::string& wellName) const
+        {
+            const auto& it = wellMap().find( wellName );
+            if (it == wellMap().end()) {
+                OPM_THROW(std::logic_error, "Could not find well " << wellName << " in well map");
+            }
+            const int well_index = it->second[0];
+            return wellIsOwned(well_index, wellName);
+        }
+
         /// The number of wells present.
         int numWells() const
         {

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -1351,6 +1351,26 @@ namespace Opm
             do_glift_optimization_ = true;
         }
 
+        int wellNameToGlobalIdx(const std::string &name) {
+            auto it = wellNameToGlobalIdx_.find(name);
+            if (it == wellNameToGlobalIdx_.end())
+                OPM_THROW(std::logic_error, "Could not find global well idx for well " << name);
+
+            return it->second;
+        }
+
+        std::string globalIdxToWellName(const int index) {
+            using Pair = std::pair<std::string, int>;
+            auto it = find_if(wellNameToGlobalIdx_.begin(), wellNameToGlobalIdx_.end(), [index](const Pair & p) {
+                return p.second == index;
+            });
+            if (it == wellNameToGlobalIdx_.end() )
+            {
+                OPM_THROW(std::logic_error, "Could not find well name for global idx " << index);
+            }
+            return it->first;
+        }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<bool> is_producer_; // Size equal to number of local wells.