Merge pull request #2709 from bska/report-guiderates

Report Guiderates at Well and Group Levels

opm/simulators/wells/BlackoilWellModel.hpp

@@ -31,6 +31,7 @@
 #include <opm/common/utility/platform_dependent/reenable_warnings.h>
 
 #include <cassert>
+#include <unordered_map>
 #include <tuple>
 
 #include <opm/parser/eclipse/EclipseState/Runspec.hpp>
@@ -195,11 +196,15 @@ namespace Opm {
             {
                 auto gvalues = ::Opm::data::GroupValues{};
 
+                const auto groupGuideRates =
+                    calculateAllGroupGuiderates(reportStepIdx, sched);
+
                 for (const auto& gname : sched.groupNames(reportStepIdx)) {
                     const auto& grup = sched.getGroup(gname, reportStepIdx);
 
                     auto& gdata = gvalues[gname];
                     this->assignGroupControl(grup, gdata);
+                    this->assignGroupGuideRates(grup, groupGuideRates, gdata);
                 }
 
                 return gvalues;
@@ -216,6 +221,9 @@ namespace Opm {
                     }
 
                     xwPos->second.current_control.isProducer = well.isProducer();
+
+                    auto& grval = xwPos->second.guide_rates;  grval.clear();
+                    grval += this->getGuideRateValues(well);
                 }
 
                 return wsrpt;
@@ -435,7 +443,20 @@ namespace Opm {
 
             void setWsolvent(const Group& group, const Schedule& schedule, const int reportStepIdx, double wsolvent);
 
+            std::unordered_map<std::string, data::GroupGuideRates>
+            calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const;
+
             void assignGroupControl(const Group& group, data::GroupData& gdata) const;
+            data::GuideRateValue getGuideRateValues(const Well& well) const;
+            data::GuideRateValue getGuideRateValues(const Group& group) const;
+            void getGuideRateValues(const GuideRate::RateVector& qs,
+                                    const bool                   is_inj,
+                                    const std::string&           wgname,
+                                    data::GuideRateValue&        grval) const;
+
+            void assignGroupGuideRates(const Group& group,
+                                       const std::unordered_map<std::string, data::GroupGuideRates>& groupGuideRates,
+                                       data::GroupData& gdata) const;
         };
 
 

opm/simulators/wells/BlackoilWellModel_impl.hpp

@@ -24,6 +24,8 @@
 #include <opm/simulators/wells/SimFIBODetails.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
 
+#include <utility>
+
 namespace Opm {
     template<typename TypeTag>
     BlackoilWellModel<TypeTag>::
@@ -2430,6 +2432,79 @@ namespace Opm {
         }
     }
 
+    template<typename TypeTag>
+    std::unordered_map<std::string, data::GroupGuideRates>
+    BlackoilWellModel<TypeTag>::
+    calculateAllGroupGuiderates(const int reportStepIdx, const Schedule& sched) const
+    {
+        auto gr = std::unordered_map<std::string, data::GroupGuideRates>{};
+        auto up = std::vector<std::string>{};
+
+        // Start at well level, accumulate contributions towards root of
+        // group tree (FIELD group).
+
+        for (const auto& wname : sched.wellNames(reportStepIdx)) {
+            if (! (this->well_state_.hasWellRates(wname) &&
+                   this->guideRate_->has(wname)))
+            {
+                continue;
+            }
+
+            const auto& well   = sched.getWell(wname, reportStepIdx);
+            const auto& parent = well.groupName();
+
+            if (parent == "FIELD") {
+                // Well parented directly to "FIELD".  Inadvisable and
+                // unexpected, but nothing to do about that here.  Just skip
+                // this guide rate contribution.
+                continue;
+            }
+
+            auto& grval = well.isInjector()
+                ? gr[parent].injection
+                : gr[parent].production;
+
+            grval += this->getGuideRateValues(well);
+            up.push_back(parent);
+        }
+
+        // Propagate accumulated guide rates up towards root of group tree.
+        // Override accumulation if there is a GUIDERAT specification that
+        // applies to a group.
+        std::sort(up.begin(), up.end());
+        auto start = 0*up.size();
+        auto u     = std::unique(up.begin(), up.end());
+        auto nu    = std::distance(up.begin(), u);
+        while (nu > 0) {
+            const auto ntot = up.size();
+
+            for (auto gi = 0*nu; gi < nu; ++gi) {
+                const auto& gname = up[start + gi];
+                const auto& group = sched.getGroup(gname, reportStepIdx);
+
+                if (this->guideRate_->has(gname)) {
+                    gr[gname].production = this->getGuideRateValues(group);
+                }
+
+                const auto parent = group.parent();
+                if (parent == "FIELD") { continue; }
+
+                gr[parent].injection  += gr[gname].injection;
+                gr[parent].production += gr[gname].production;
+                up.push_back(parent);
+            }
+
+            start = ntot;
+
+            auto begin = up.begin() + ntot;
+            std::sort(begin, up.end());
+            u  = std::unique(begin, up.end());
+            nu = std::distance(begin, u);
+        }
+
+        return gr;
+    }
+
     template<typename TypeTag>
     void
     BlackoilWellModel<TypeTag>::
@@ -2466,4 +2541,114 @@ namespace Opm {
         }
     }
 
+    template <typename TypeTag>
+    data::GuideRateValue
+    BlackoilWellModel<TypeTag>::
+    getGuideRateValues(const Well& well) const
+    {
+        auto grval = data::GuideRateValue{};
+
+        assert (this->guideRate_ != nullptr);
+
+        const auto& wname = well.name();
+        if (! this->well_state_.hasWellRates(wname)) {
+            // No flow rates for 'wname' -- might be before well comes
+            // online (e.g., for the initial condition before simulation
+            // starts).
+            return grval;
+        }
+
+        if (! this->guideRate_->has(wname)) {
+            // No guiderates exist for 'wname'.
+            return grval;
+        }
+
+        const auto qs = WellGroupHelpers::
+            getRateVector(this->well_state_, this->phase_usage_, wname);
+
+        this->getGuideRateValues(qs, well.isInjector(), wname, grval);
+
+        return grval;
+    }
+
+    template <typename TypeTag>
+    data::GuideRateValue
+    BlackoilWellModel<TypeTag>::
+    getGuideRateValues(const Group& group) const
+    {
+        auto grval = data::GuideRateValue{};
+
+        assert (this->guideRate_ != nullptr);
+
+        const auto& gname = group.name();
+        if (! this->well_state_.hasProductionGroupRates(gname)) {
+            // No flow rates for 'gname' -- might be before group comes
+            // online (e.g., for the initial condition before simulation
+            // starts).
+            return grval;
+        }
+
+        if (! this->guideRate_->has(gname)) {
+            // No guiderates exist for 'gname'.
+            return grval;
+        }
+
+        const auto qs = WellGroupHelpers::
+            getProductionGroupRateVector(this->well_state_, this->phase_usage_, gname);
+
+        const auto is_inj = false; // This procedure only applies to G*PGR.
+        this->getGuideRateValues(qs, is_inj, gname, grval);
+
+        return grval;
+    }
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    getGuideRateValues(const GuideRate::RateVector& qs,
+                       const bool                   is_inj,
+                       const std::string&           wgname,
+                       data::GuideRateValue&        grval) const
+    {
+        auto getGR = [this, &wgname, &qs](const GuideRateModel::Target t)
+        {
+            return this->guideRate_->get(wgname, t, qs);
+        };
+
+        // Note: GuideRate does currently (2020-07-20) not support Target::RES.
+        grval.set(data::GuideRateValue::Item::Gas,
+                  getGR(GuideRateModel::Target::GAS));
+
+        grval.set(data::GuideRateValue::Item::Water,
+                  getGR(GuideRateModel::Target::WAT));
+
+        if (! is_inj) {
+            // Producer.  Extract "all" guiderate values.
+            grval.set(data::GuideRateValue::Item::Oil,
+                      getGR(GuideRateModel::Target::OIL));
+        }
+    }
+
+    template <typename TypeTag>
+    void
+    BlackoilWellModel<TypeTag>::
+    assignGroupGuideRates(const Group& group,
+                          const std::unordered_map<std::string, data::GroupGuideRates>& groupGuideRates,
+                          data::GroupData& gdata) const
+    {
+        auto& prod = gdata.guideRates.production;  prod.clear();
+        auto& inj  = gdata.guideRates.injection;   inj .clear();
+
+        auto xgrPos = groupGuideRates.find(group.name());
+        if ((xgrPos == groupGuideRates.end()) ||
+            !this->guideRate_->has(group.name()))
+        {
+            // No guiderates defined for this group.
+            return;
+        }
+
+        const auto& xgr = xgrPos->second;
+        prod = xgr.production;
+        inj  = xgr.injection;
+    }
 } // namespace Opm

opm/simulators/wells/WellGroupHelpers.cpp

@@ -24,6 +24,28 @@
 #include <algorithm>
 #include <vector>
 
+namespace {
+    Opm::GuideRate::RateVector
+    getGuideRateVector(const std::vector<double>& rates, const Opm::PhaseUsage& pu)
+    {
+        using Opm::BlackoilPhases;
+
+        double oilRate = 0.0;
+        if (pu.phase_used[BlackoilPhases::Liquid])
+            oilRate = rates[pu.phase_pos[BlackoilPhases::Liquid]];
+
+        double gasRate = 0.0;
+        if (pu.phase_used[BlackoilPhases::Vapour])
+            gasRate = rates[pu.phase_pos[BlackoilPhases::Vapour]];
+
+        double waterRate = 0.0;
+        if (pu.phase_used[BlackoilPhases::Aqua])
+            waterRate = rates[pu.phase_pos[BlackoilPhases::Aqua]];
+
+        return {oilRate, gasRate, waterRate};
+    }
+} // namespace Anonymous
+
 namespace Opm
 {
 
@@ -566,22 +588,14 @@ namespace WellGroupHelpers
     GuideRate::RateVector
     getRateVector(const WellStateFullyImplicitBlackoil& well_state, const PhaseUsage& pu, const std::string& name)
     {
-        const std::vector<double>& rates = well_state.currentWellRates(name);
-        double oilRate = 0.0;
-        if (pu.phase_used[BlackoilPhases::Liquid])
-            oilRate = rates[pu.phase_pos[BlackoilPhases::Liquid]];
-
-        double gasRate = 0.0;
-        if (pu.phase_used[BlackoilPhases::Vapour])
-            gasRate = rates[pu.phase_pos[BlackoilPhases::Vapour]];
-
-        double waterRate = 0.0;
-        if (pu.phase_used[BlackoilPhases::Aqua])
-            waterRate = rates[pu.phase_pos[BlackoilPhases::Aqua]];
-
-        return GuideRate::RateVector {oilRate, gasRate, waterRate};
+        return getGuideRateVector(well_state.currentWellRates(name), pu);
     }
 
+    GuideRate::RateVector
+    getProductionGroupRateVector(const WellStateFullyImplicitBlackoil& well_state, const PhaseUsage& pu, const std::string& group_name)
+    {
+        return getGuideRateVector(well_state.currentProductionGroupRates(group_name), pu);
+    }
 
     double getGuideRate(const std::string& name,
                         const Schedule& schedule,
@@ -699,7 +713,6 @@ namespace WellGroupHelpers
         return num_wells;
     }
 
-
     FractionCalculator::FractionCalculator(const Schedule& schedule,
                                            const WellStateFullyImplicitBlackoil& well_state,
                                            const int report_step,
@@ -791,21 +804,7 @@ namespace WellGroupHelpers
 
     GuideRate::RateVector FractionCalculator::getGroupRateVector(const std::string& group_name)
     {
-
-        std::vector<double> groupRates = well_state_.currentProductionGroupRates(group_name);
-        double oilRate = 0.0;
-        if (pu_.phase_used[BlackoilPhases::Liquid])
-            oilRate = groupRates[pu_.phase_pos[BlackoilPhases::Liquid]];
-
-        double gasRate = 0.0;
-        if (pu_.phase_used[BlackoilPhases::Vapour])
-            gasRate = groupRates[pu_.phase_pos[BlackoilPhases::Vapour]];
-
-        double waterRate = 0.0;
-        if (pu_.phase_used[BlackoilPhases::Aqua])
-            waterRate = groupRates[pu_.phase_pos[BlackoilPhases::Aqua]];
-
-        return GuideRate::RateVector {oilRate, gasRate, waterRate};
+        return getProductionGroupRateVector(this->well_state_, this->pu_, group_name);
     }
 
 

opm/simulators/wells/WellGroupHelpers.hpp

@@ -261,6 +261,8 @@ namespace WellGroupHelpers
     GuideRate::RateVector
     getRateVector(const WellStateFullyImplicitBlackoil& well_state, const PhaseUsage& pu, const std::string& name);
 
+    GuideRate::RateVector
+    getProductionGroupRateVector(const WellStateFullyImplicitBlackoil& well_state, const PhaseUsage& pu, const std::string& group_name);
 
     double getGuideRate(const std::string& name,
                         const Schedule& schedule,
@@ -280,7 +282,6 @@ namespace WellGroupHelpers
                            const Phase& injectionPhase,
                            const PhaseUsage& pu);
 
-
     int groupControlledWells(const Schedule& schedule,
                              const WellStateFullyImplicitBlackoil& well_state,
                              const int report_step,
@@ -314,7 +315,6 @@ namespace WellGroupHelpers
         PhaseUsage pu_;
     };
 
-    GuideRate::RateVector getGroupRateVector(const std::string& group_name);
 
 
     double fractionFromGuideRates(const std::string& name,

opm/simulators/wells/WellStateFullyImplicitBlackoil.hpp

@@ -379,6 +379,10 @@ namespace Opm
             return it->second;
         }
 
+        bool hasWellRates(const std::string& wellName) const {
+            return this->well_rates.find(wellName) != this->well_rates.end();
+        }
+
         void setCurrentProductionGroupRates(const std::string& groupName, const std::vector<double>& rates ) {
             production_group_rates[groupName] = rates;
         }
@@ -392,6 +396,10 @@ namespace Opm
             return it->second;
         }
 
+        bool hasProductionGroupRates(const std::string& groupName) const {
+            return this->production_group_rates.find(groupName) != this->production_group_rates.end();
+        }
+        
         void setCurrentProductionGroupReductionRates(const std::string& groupName, const std::vector<double>& target ) {
             production_group_reduction_rates[groupName] = target;
         }
@@ -1066,9 +1074,9 @@ namespace Opm
             auto it = wellNameToGlobalIdx_.find(name);
 
             if (it == wellNameToGlobalIdx_.end())
-                OPM_THROW(std::logic_error, "Could not find global injection group for well" << name);
+                OPM_THROW(std::logic_error, "Could not find global injection group for well " << name);
 
-            return globalIsInjectionGrup_[it->second];
+            return globalIsInjectionGrup_[it->second] != 0;
         }
 
         bool isProductionGrup(const std::string& name) const {
@@ -1076,9 +1084,9 @@ namespace Opm
             auto it = wellNameToGlobalIdx_.find(name);
 
             if (it == wellNameToGlobalIdx_.end())
-                OPM_THROW(std::logic_error, "Could not find global injection group for well" << name);
+                OPM_THROW(std::logic_error, "Could not find global production group for well " << name);
 
-            return globalIsProductionGrup_[it->second];
+            return globalIsProductionGrup_[it->second] != 0;
         }
 
     private: