Completly revamped the way group control is checked

opm/core/WellCollection.cpp

@@ -93,15 +93,17 @@ namespace Opm
         return NULL;
     }
     
-    void WellCollection::conditionsMet(const std::vector<double>& well_bhp,
+    bool WellCollection::conditionsMet(const std::vector<double>& well_bhp,
                                        const std::vector<double>& well_rate, 
-                                       const UnstructuredGrid& grid,
+                                       double epsilon)
-                                       WellControlResult& result,
-                                       double epsilon) const
     {
-        for (size_t i = 0; i < leaf_nodes_.size(); i++) {
+        for (size_t i = 0; i < roots_.size(); i++) {
-            leaf_nodes_[i]->conditionsMet(well_bhp, well_rate, grid, result, epsilon);
+            WellPhasesSummed phases;
+            if(!roots_[i]->conditionsMet(well_bhp, well_rate, phases, epsilon)) {
+                return false;
+            }
         }
+        return true;
     }
 
     void WellCollection::calculateGuideRates()
@@ -111,7 +113,7 @@ namespace Opm
         }
     }
     
-    void WellCollection::setWellsPointer(const Wells* wells) {
+    void WellCollection::setWellsPointer(Wells* wells) {
         for(size_t i = 0; i < leaf_nodes_.size(); i++) {
             leaf_nodes_[i]->setWellsPointer(wells, i);
         }

opm/core/WellCollection.hpp

@@ -43,14 +43,12 @@ namespace Opm
                       const std::string& parent,
                       const EclipseGridParser& deck);
 
-        /// Builds the WellControlResult object for the current well group hierachy. 
+        bool conditionsMet(const std::vector<double>& well_bhp,
-        void conditionsMet(const std::vector<double>& well_bhp,
                            const std::vector<double>& well_rate, 
-                           const UnstructuredGrid& grid,
+                           const double epsilon=1e-8);
-                           WellControlResult& result, 
+        
-                           const double epsilon=1e-8) const;
         /// Adds the well pointer to each leaf node (does not take ownership).
-        void setWellsPointer(const Wells* wells);
+        void setWellsPointer(Wells* wells);
         
         const std::vector<WellNode*>& getLeafNodes() const;
         

opm/core/WellsGroup.cpp

@@ -12,6 +12,18 @@
 namespace Opm
 {
 
+    WellPhasesSummed::WellPhasesSummed() 
+    : bhp_sum(0.0), rate_sum(0.0)
+    {
+        
+    }
+    
+    void WellPhasesSummed::operator+=(const WellPhasesSummed& other) 
+    {
+        rate_sum += other.rate_sum;
+        bhp_sum += other.bhp_sum;
+    }
+    
     WellsGroupInterface::WellsGroupInterface(const std::string& myname,
             ProductionSpecification prod_spec,
             InjectionSpecification inje_spec)
@@ -115,64 +127,79 @@ namespace Opm
         }
     }
     
+    void WellsGroup::applyControl(const WellControlType type) 
+    {
+        for (size_t i = 0; i < children_.size(); ++i) {
+            children_[i]->applyControl(type);
+        }
+    }
 
-    void WellsGroup::conditionsMet(const std::vector<double>& well_bhp,
+    
+    bool WellsGroup::conditionsMet(const std::vector<double>& well_bhp, 
                                    const std::vector<double>& well_rate, 
-                                   const UnstructuredGrid& grid,
+                                   WellPhasesSummed& summed_phases, 
-                                   const struct Wells* wells,
-                                   int index_of_well,
-                                   WellControlResult& result,
                                    const double epsilon)
     {
-        if (parent_ != NULL) {
+        WellPhasesSummed child_phases_summed;
-            (static_cast<WellsGroup*> (parent_))->conditionsMet(well_bhp, 
+        for(size_t i = 0; i < children_.size(); ++i) {
-                    well_rate,grid, wells, index_of_well, result, epsilon);
+            WellPhasesSummed current_child_phases_summed;
+            if(!children_[i]->conditionsMet(well_bhp, well_rate, 
+                                            current_child_phases_summed, epsilon)) {
+                return false;
+            }
+            child_phases_summed += current_child_phases_summed;
         }
         
-        int number_of_leaf_nodes = numberOfLeafNodes();
 
-        double bhp_target = 1e100;
+        double bhp_target = std::min(injSpec().BHP_limit_, prodSpec().BHP_limit_);
-        double rate_target = 1e100;
+        double rate_target = std::min(injSpec().fluid_volume_max_rate_, 
-        switch(wells->type[index_of_well]) {
+                                      prodSpec().fluid_volume_max_rate_);
-        case INJECTOR:
-        {
-            const InjectionSpecification& inje_spec = injSpec();
-            bhp_target = inje_spec.BHP_limit_ / number_of_leaf_nodes;
-            rate_target = inje_spec.fluid_volume_max_rate_ / number_of_leaf_nodes;
-            break;
-        }
-        case PRODUCER:
-        {
-            const ProductionSpecification& prod_spec = prodSpec();
-            bhp_target = prod_spec.BHP_limit_ / number_of_leaf_nodes;
-            rate_target = prod_spec.fluid_volume_max_rate_ / number_of_leaf_nodes;
-            break;
-        }
-        }
         
-        if (well_bhp[index_of_well] - bhp_target > epsilon) {
+        double bhp_sum = child_phases_summed.bhp_sum;
+        double rate_sum = child_phases_summed.rate_sum;
+        if (bhp_sum - bhp_target > epsilon) {
             std::cout << "BHP not met" << std::endl;
             std::cout << "BHP limit was " << bhp_target << std::endl;
-            std::cout << "Actual bhp was " << well_bhp[index_of_well] << std::endl;
+            std::cout << "Actual bhp was " << bhp_sum << std::endl;
-            ExceedInformation info;
-            info.group_name_ = name();
-            info.surplus_ = well_bhp[index_of_well] - bhp_target;
-            info.well_index_ = index_of_well;
-            result.bhp_.push_back(info);
             
+            switch(prodSpec().procedure_) {
+            case ProductionSpecification::WELL:
+                getWorstOffending(well_bhp).first->shutWell();
+                return false;
+                break;
+            case ProductionSpecification::RATE:
+                applyControl(BHP);
+                return false;
+                break;
+            default:
+                // Nothing do to;
+                break;
+            }
         }
-        if(well_rate[index_of_well] - rate_target > epsilon) {
+        if(rate_sum - rate_target > epsilon) {
             std::cout << "well_rate not met" << std::endl;
-            std::cout << "target = " << rate_target << ", well_rate[index_of_well] = " << well_rate[index_of_well] << std::endl;
+            std::cout << "target = " << rate_target 
+                      << ", well_rate[index_of_well] = " 
+                      << rate_sum << std::endl;
             std::cout << "Group name = " << name() << std::endl;
             
-            ExceedInformation info;
+            switch(prodSpec().procedure_) {
-            info.group_name_ = name();
+            case ProductionSpecification::WELL:
-            info.surplus_ = well_rate[index_of_well] - rate_target;
+                getWorstOffending(well_rate).first->shutWell();
-            info.well_index_ = index_of_well;
+                return false;
-            result.fluid_rate_.push_back(info);
+                break;
+            case ProductionSpecification::RATE:
+                applyControl(RATE);
+                return false;
+                break;
+            default:
+                // Nothing do to;
+                break;
+            }
         }
         
+        summed_phases += child_phases_summed;
+        return true;
     }
 
     void WellsGroup::addChild(std::tr1::shared_ptr<WellsGroupInterface> child)
@@ -193,6 +220,18 @@ namespace Opm
         return sum;
     }
     
+    std::pair<WellNode*, double> WellsGroup::getWorstOffending(const std::vector<double>& values) 
+    {
+        std::pair<WellNode*, double> max;
+        for (size_t i = 0; i < children_.size(); i++) {
+            std::pair<WellNode*, double> child_max = children_[i]->getWorstOffending(values);
+            if (i == 0 || max.second < child_max.second) {
+                max = child_max;
+            }
+        }
+        return max;
+    }
+    
     WellNode::WellNode(const std::string& myname,
             ProductionSpecification prod_spec,
             InjectionSpecification inj_spec)
@@ -200,18 +239,12 @@ namespace Opm
     {
     }
 
-    void WellNode::conditionsMet(const std::vector<double>& well_bhp,
+    bool WellNode::conditionsMet(const std::vector<double>& well_bhp, 
                                  const std::vector<double>& well_rate, 
-                                 const UnstructuredGrid& grid,
+                                 WellPhasesSummed& summed_phases, 
-                                 WellControlResult& result,
                                  const double epsilon)
     {
         
-        if (parent_ != NULL) {
-            (static_cast<WellsGroup*> (parent_))
-                ->conditionsMet(well_bhp, well_rate, grid, wells_,
-                                self_index_, result, epsilon);
-        }
         
         // Check for self:
         if (wells_->type[self_index_] == PRODUCER) {
@@ -223,21 +256,14 @@ namespace Opm
                 std::cout << "BHP exceeded, bhp_diff = " << bhp_diff << std::endl;
                 std::cout << "BHP_limit = " << prodSpec().BHP_limit_ << std::endl;
                 std::cout << "BHP = " << well_bhp[self_index_] << std::endl; 
-
+                shutWell();
-                ExceedInformation info;
+                return false;
-                info.group_name_ = name();
-                info.surplus_ = bhp_diff;
-                info.well_index_ = self_index_;
-                result.bhp_.push_back(info);
             }
             
             if (rate_diff > epsilon) {
-                ExceedInformation info;
-                info.group_name_ = name();
-                info.surplus_ = rate_diff;
-                info.well_index_ = self_index_;
-                result.fluid_rate_.push_back(info);
                 std::cout << "Rate exceeded, rate_diff = " << rate_diff << std::endl;
+                shutWell();
+                return false;
             }
         } else {
             double bhp_diff = well_bhp[self_index_] - injSpec().BHP_limit_;
@@ -246,22 +272,20 @@ namespace Opm
             
            if (bhp_diff > epsilon) {
                std::cout << "BHP exceeded, bhp_diff = " << bhp_diff<<std::endl;
-               ExceedInformation info;
+               shutWell();
-               info.group_name_ = name();
+               return false;
-               info.surplus_ = bhp_diff;
-               info.well_index_ = self_index_;
-               result.bhp_.push_back(info);
            }
            if (rate_diff > epsilon) {
                std::cout << "Flow diff exceeded, flow_diff = " << rate_diff << std::endl;
-               ExceedInformation info;
+               shutWell();
-               info.group_name_ = name();
+               return false;
-               info.surplus_ = rate_diff;
-               info.well_index_ = self_index_;
-               result.fluid_rate_.push_back(info);
            }
             
         }
+        
+        summed_phases.bhp_sum = well_bhp[self_index_];
+        summed_phases.rate_sum = well_rate[self_index_];
+        return true;
     }
 
     WellsGroupInterface* WellNode::findGroup(std::string name_of_node)
@@ -279,7 +303,7 @@ namespace Opm
         return true;
     }
 
-    void WellNode::setWellsPointer(const struct Wells* wells, int self_index)
+    void WellNode::setWellsPointer(Wells* wells, int self_index)
     {
         wells_ = wells;
         self_index_ = self_index;
@@ -295,6 +319,40 @@ namespace Opm
         return 1;
     }
     
+    void WellNode::shutWell() 
+    {
+        wells_->ctrls[self_index_]->target[0] = 0.0;
+    }
+    
+    std::pair<WellNode*, double> WellNode::getWorstOffending(const std::vector<double>& values) {
+        return std::make_pair<WellNode*, double>(this, values[self_index_]);
+    }
+    
+    void WellNode::applyControl(const WellControlType type)
+    {
+        wells_->ctrls[self_index_]->type[0] = type;
+        double target = 0.0;
+        switch(type) {
+        case BHP:
+            if(wells_->type[self_index_] == INJECTOR) {
+                target = injSpec().BHP_limit_;
+            }
+            else {
+                target = prodSpec().BHP_limit_;
+            }
+            break;
+        case RATE:
+            if(wells_->type[self_index_] == INJECTOR) {
+                target = injSpec().fluid_volume_max_rate_;
+            }
+            else {
+                target = prodSpec().fluid_volume_max_rate_;
+            }
+            break;
+        }
+        wells_->ctrls[self_index_]->target[0] = target;
+    }
+
     namespace
     {
 

opm/core/WellsGroup.hpp

@@ -10,17 +10,13 @@
 
 namespace Opm
 {
+    class WellNode;
+    struct WellPhasesSummed {
+        WellPhasesSummed();
+        double bhp_sum;
+        double rate_sum;
         
-    struct ExceedInformation {
+        void operator+=(const WellPhasesSummed& other);
-        std::string group_name_;
-        int well_index_;
-        double surplus_;
-    };
-    
-    struct WellControlResult {
-        std::vector<ExceedInformation> oil_rate_;
-        std::vector<ExceedInformation> fluid_rate_;
-        std::vector<ExceedInformation> bhp_;
     };
 
     class WellsGroupInterface
@@ -64,6 +60,17 @@ namespace Opm
         /// Calculates the number of leaf nodes in the given group. 
         /// A leaf node is defined to have one leaf node in its group.
         virtual int numberOfLeafNodes() = 0;
+
+        /// Fills the WellControlResult parameter with all exceed information
+        virtual bool conditionsMet(const std::vector<double>& well_bhp,
+                                   const std::vector<double>& well_rate,
+                                   WellPhasesSummed& summed_phases,
+                                   const double epsilon = 1e-8) = 0;
+        
+        virtual void applyControl(const WellControlType type) = 0;
+        
+        virtual std::pair<WellNode*, double> getWorstOffending(const std::vector<double>& values) = 0;
+        
     protected:
            WellsGroupInterface* parent_;
 
@@ -86,17 +93,18 @@ namespace Opm
 
         void addChild(std::tr1::shared_ptr<WellsGroupInterface> child);
         
-        void conditionsMet(const std::vector<double>& well_bhp,
+        virtual bool conditionsMet(const std::vector<double>& well_bhp,
-                           const std::vector<double>& well_rate,
+                                   const std::vector<double>& well_rate,
-                           const UnstructuredGrid& grid,
+                                   WellPhasesSummed& summed_phases,
-                           const struct Wells* wells,
+                                   const double epsilon = 1e-8);
-                           int index_of_well,
+        
-                           WellControlResult& result,
-                           double epsilon = 1e-8);
         
         virtual void calculateGuideRates();
 
         virtual int numberOfLeafNodes();
+        virtual std::pair<WellNode*, double> getWorstOffending(const std::vector<double>& values);
+        virtual void applyControl(const WellControlType type);
+
     private:
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > children_;
     };
@@ -111,19 +119,25 @@ namespace Opm
                 InjectionSpecification inj_spec);
 
         virtual WellsGroupInterface* findGroup(std::string name_of_node);
-        virtual void conditionsMet(const std::vector<double>& well_bhp,
+        virtual bool conditionsMet(const std::vector<double>& well_bhp,
                                    const std::vector<double>& well_rate,
-                                   const UnstructuredGrid& grid,
+                                   WellPhasesSummed& summed_phases,
-                                   WellControlResult& result,
+                                   const double epsilon = 1e-8);
-                                   double epsilon=1e-8);
+        
         virtual bool isLeafNode() const;
         
-        void setWellsPointer(const struct Wells* wells, int self_index);
+        void setWellsPointer(Wells* wells, int self_index);
         
         virtual void calculateGuideRates();
         virtual int numberOfLeafNodes();
+        
+        void shutWell();
+        
+        virtual std::pair<WellNode*, double> getWorstOffending(const std::vector<double>& values);
+        virtual void applyControl(const WellControlType type);
+        
     private:
-        const struct Wells* wells_;
+        Wells* wells_;
         int self_index_;
     };
 

opm/core/WellsManager.cpp

@@ -637,120 +637,10 @@ namespace Opm
         return well_collection_;
     }
 
-    
+    bool WellsManager::conditionsMet(const std::vector<double>& well_bhp, 
-    /// Apply control results
+                                     const std::vector<double>& well_rate) 
-    /// \param[in] result The result of a run to conditionsMet on WellCollection
-    /// \param[in] wellCollection The WellCollection on which the control is to be issued on
-    void WellsManager::applyControl(const WellControlResult& result) 
     {
-        // Check oil
+        return well_collection_.conditionsMet(well_bhp, well_rate);
-        std::map<std::string, std::vector<ExceedInformation> > oil_exceed;
-        for(size_t i = 0; i < result.oil_rate_.size(); i++) {
-            oil_exceed[result.oil_rate_[i].group_name_].push_back(result.oil_rate_[i]);
-        }
-        
-        applyControl(oil_exceed, ProductionSpecification::ORAT);
-        
-        
-         // Check fluid
-        std::map<std::string, std::vector<ExceedInformation> > fluid_exceed;
-        for(size_t i = 0; i < result.fluid_rate_.size(); i++) {
-            fluid_exceed[result.oil_rate_[i].group_name_].push_back(result.fluid_rate_[i]);
-        }
-        
-        applyControl(fluid_exceed, ProductionSpecification::LRAT);
-        
-        // Check BHP
-        std::map<std::string, std::vector<ExceedInformation> > bhp_exceed;
-        for(size_t i = 0; i < result.bhp_.size(); i++) {
-            bhp_exceed[result.oil_rate_[i].group_name_].push_back(result.bhp_[i]);
-        }
-        
-        applyControl(fluid_exceed, ProductionSpecification::BHP); 
-        
-        
-        // Apply guide rates:
-        for (int i = 0; i < w_->number_of_wells; i++) {
-            if (well_collection_.getLeafNodes()[i]->prodSpec().control_mode_ == ProductionSpecification::GRUP) {
-                switch (well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_type_) {
-                case ProductionSpecification::OIL:
-                {
-                    // Not handled at the moment
-                }
-                case ProductionSpecification::NONE_GRT:
-                {
-                    // Will use the group control type:
-                    const ProductionSpecification& parent_prod_spec =
-                            well_collection_.getLeafNodes()[i]->getParent()->prodSpec();
-                    double guide_rate = well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_;
-                    if (parent_prod_spec.control_mode_ == ProductionSpecification::LRAT) {
-                        w_->ctrls[i]->target[0] = guide_rate * parent_prod_spec.liquid_max_rate_;
-                        w_->ctrls[i]->type[0] = RATE;
-                    } else {
-                        THROW("Unhandled group control mode " << parent_prod_spec.control_mode_);
-                    }
-
-                }
-                default:
-                    // Do nothing
-                    break;
-                }
-            }
-        }
-    }
-
-    /// Apply control results for a specific target (OIL, WATER, etc)
-    /// \param[in] exceed_info will for each group name contain all the 
-    ///                        exceed informations for the given mode.
-    /// \param[in] well_collection The associated well_collection.
-    /// \param[in] mode The ControlMode to which the violations apply.
-    void WellsManager::applyControl(const std::map<std::string, std::vector<ExceedInformation> >& exceed_info,
-                                    ProductionSpecification::ControlMode mode) 
-    {
-        std::map<std::string, std::vector<ExceedInformation> >::const_iterator it;
-        
-        for(it = exceed_info.begin(); it != exceed_info.end(); ++it) {
-
-            
-            std::string group_name = it->first;
-            
-            WellsGroupInterface* group = well_collection_.findNode(group_name);
-            if(group->isLeafNode()) {
-                // Just shut the well
-                int well_index = it->second[0].well_index_;
-                w_->ctrls[well_index]->target[0] = 0.0;
-            }
-            else {
-                switch(group->prodSpec().procedure_) {
-                case ProductionSpecification::WELL:
-                {
-                    // Shut the worst offending well
-                    double max_exceed = 0.0;
-                    int exceed_index = -1;
-                    for(size_t i = 0; i < it->second.size(); i++) {
-                        if(max_exceed <= it->second[i].surplus_) {
-                            exceed_index = it->second[i].well_index_;
-                            max_exceed = it->second[i].surplus_;
-                        }
-                    }
-                    
-                    w_->ctrls[exceed_index]->target[0] = 0.0;
-                    break;
-                }
-                
-                case ProductionSpecification::RATE:
-                {
-                    // Now we need to set the group control mode to the active one
-                    group->prodSpec().control_mode_ = mode;
-                    break;
-                }
-                
-                default:
-                    // Do nothing for now
-                    break;
-                }
-            }
-        }
     }
 
 } // namespace Opm

opm/core/WellsManager.hpp

@@ -51,9 +51,6 @@ namespace Opm
 	/// Destructor.
 	~WellsManager();
         
-        /// Apply control results
-        /// \param[in] result The result of a run to conditionsMet on WellCollection
-        void applyControl(const WellControlResult& result);
         
 	/// Access the managed Wells.
 	/// The method is named similarly to c_str() in std::string,
@@ -63,6 +60,9 @@ namespace Opm
         /// Access the well group hierarchy.
         const WellCollection& wellCollection() const;
         
+        bool conditionsMet(const std::vector<double>& well_bhp, 
+                           const std::vector<double>& well_rate);
+
     private:
 	// Disable copying and assignment.
 	WellsManager(const WellsManager& other);
@@ -72,13 +72,7 @@ namespace Opm
 	Wells* w_;
         WellCollection well_collection_;
 
-        /// Apply control results for a specific target (OIL, WATER, etc)
+       
-        /// \param[in] exceed_info will for each group name contain all the 
-        ///                        exceed informations for the given mode.
-        /// \param[in] well_collection The associated well_collection.
-        /// \param[in] mode The ControlMode to which the violations apply.
-        void applyControl(const std::map<std::string, std::vector<ExceedInformation> >& exceed_info,
-                          ProductionSpecification::ControlMode mode);
 
 
     };