Added a (small) hack to wellmanager (defaulting injected_phase to water)

opm/core/ProductionSpecification.hpp

@@ -11,7 +11,7 @@ namespace Opm
 
         enum ControlMode
         {
-            NONE_CM, ORAT, WRAT, LRAT, REIN, RESV, VREP, WGRA, FLD, GRUP
+            NONE_CM, ORAT, WRAT, LRAT, REIN, RESV, VREP, WGRA, FLD, GRUP, BHP
         };
 
         enum Procedure

opm/core/WellCollection.cpp

@@ -88,17 +88,24 @@ namespace Opm
         return NULL;
     }
     
-    bool WellCollection::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
+    const WellsGroupInterface* WellCollection::findNode(std::string name) const
+    {
+
+        for (size_t i = 0; i < roots_.size(); i++) {
+            WellsGroupInterface* result = roots_[i]->findGroup(name);
+            if (result) {
+                return result;
+            }
+        }
+        return NULL;
+    }
+    
+    void WellCollection::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
                            const UnstructuredGrid& grid, const std::vector<double>& saturations, 
                            WellControlResult& result, double epsilon) const {
         for(size_t i = 0; i < leaf_nodes_.size(); i++) {
-            if(! static_cast<WellNode*>(leaf_nodes_[i].get())->conditionsMet(well_bhp, well_rate, grid, saturations, result, epsilon) ) {
+            static_cast<WellNode*>(leaf_nodes_[i].get())->conditionsMet(well_bhp, well_rate, grid, saturations, result, epsilon);
-                return false;
+        }        
-            }
-        }
-        
-        return true;
-        
     }
 
     void WellCollection::calculateGuideRates()

opm/core/WellCollection.hpp

@@ -40,7 +40,7 @@ namespace Opm
                 const EclipseGridParser& deck);
 
         
-        bool conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
+        void conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
                            const UnstructuredGrid& grid, const std::vector<double>& saturations, 
                            WellControlResult& result, 
                            double epsilon=1e-8) const;
@@ -48,6 +48,9 @@ namespace Opm
         const std::vector<std::tr1::shared_ptr<WellsGroupInterface> >& getLeafNodes() const;
         
         void calculateGuideRates();
+        
+        WellsGroupInterface* findNode(std::string name);
+        const WellsGroupInterface* findNode(std::string name) const;
     private:
         // To account for the possibility of a forest
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > roots_;
@@ -55,7 +58,7 @@ namespace Opm
         // This will be used to traverse the bottom nodes.
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > leaf_nodes_;
         
-        WellsGroupInterface* findNode(std::string name);
+        
     };
 
 } // namespace Opm

opm/core/WellsGroup.cpp

@@ -116,23 +116,18 @@ namespace Opm
     }
 
     
-    bool WellsGroup::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
+    void WellsGroup::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
                                    const UnstructuredGrid& grid, const std::vector<double>& saturations,
                                    const struct Wells* wells, int index_of_well, WellControlResult& result,
                                    double epsilon)
     {
         if (parent_ != NULL) {
-            bool parent_ok = 
+            (static_cast<WellsGroup*> (parent_))->conditionsMet(well_bhp, 
-                (static_cast<WellsGroup*> (parent_))->conditionsMet(well_bhp, 
                     well_rate,grid, saturations, wells, index_of_well, result, epsilon);
-            if (!parent_ok) {
-                return false;
-            }
         }
         
         int number_of_leaf_nodes = numberOfLeafNodes();
 
-        bool shut_down_on_exceed = false;
         double bhp_target = 1e100;
         double rate_target = 1e100;
         switch(wells->type[index_of_well]) {
@@ -148,7 +143,6 @@ namespace Opm
             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;
-            shut_down_on_exceed = prodSpec().procedure_ == ProductionSpecification::WELL;
             break;
         }
         }
@@ -156,30 +150,26 @@ namespace Opm
         if (well_bhp[index_of_well] - 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 " << well_bhp[index_of_well] << 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);
             
-            if(shut_down_on_exceed) {
-                // Shut down well
-                // Dirty hack for now
-                struct Wells* non_const_wells = const_cast<struct Wells*>(wells);
-                non_const_wells->ctrls[index_of_well]->target[0] = 0.0;
-            }
-            return false;
         }
         if(well_rate[index_of_well] - 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 << "Group name = " << name() << std::endl;
             
-            if(shut_down_on_exceed) {
+            ExceedInformation info;
-                // Shut down well
+            info.group_name_ = name();
-                // Dirty hack for now
+            info.surplus_ = well_rate[index_of_well] - rate_target;
-                struct Wells* non_const_wells = const_cast<struct Wells*>(wells);
+            info.well_index_ = index_of_well;
-                non_const_wells->ctrls[index_of_well]->target[0] = 0.0;
+            result.fluid_rate_.push_back(info);
-            }
-            return false;
         }
-        return true;
+
     }
 
     void WellsGroup::addChild(std::tr1::shared_ptr<WellsGroupInterface> child)
@@ -207,23 +197,20 @@ namespace Opm
     {
     }
 
-    bool WellNode::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
+    void WellNode::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
                                  const UnstructuredGrid& grid, const std::vector<double>& saturations,
                                  WellControlResult& result, double epsilon)
     {
         
         if (parent_ != NULL) {
-            bool parent_ok = (static_cast<WellsGroup*> (parent_))->conditionsMet(well_bhp, 
+         (static_cast<WellsGroup*> (parent_))->conditionsMet(well_bhp, 
-                                                                                 well_rate, 
+                                                             well_rate, 
-                                                                                 grid, 
+                                                             grid, 
-                                                                                 saturations, 
+                                                             saturations, 
-                                                                                 wells_,
+                                                             wells_,
-                                                                                 self_index_,
+                                                             self_index_,
-                                                                                 result, 
+                                                             result, 
-                                                                                 epsilon);
+                                                             epsilon);
-            if (!parent_ok) {
-                return false;
-            }
         }
 
         // Check for self:
@@ -237,29 +224,44 @@ namespace Opm
                 std::cout << "BHP_limit = " << prodSpec().BHP_limit_ << std::endl;
                 std::cout << "BHP = " << well_bhp[self_index_] << std::endl; 
 
-                return false;
+                ExceedInformation info;
+                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;
-                return false;
             }
         } else {
             double bhp_diff = well_bhp[self_index_] - injSpec().BHP_limit_;
-            double flow_diff = well_rate[self_index_] - injSpec().fluid_volume_max_rate_;
+            double rate_diff = well_rate[self_index_] - injSpec().fluid_volume_max_rate_;
             
             
            if(bhp_diff > epsilon) {
                std::cout << "BHP exceeded, bhp_diff = " << bhp_diff<<std::endl;
-               return false;
+               ExceedInformation info;
+               info.group_name_ = name();
+               info.surplus_ = bhp_diff;
+               info.well_index_ = self_index_;
+               result.bhp_.push_back(info);
            }
-            if(flow_diff > epsilon) {
+            if(rate_diff > epsilon) {
-                std::cout << "Flow diff exceeded, flow_diff = " << flow_diff << std::endl;
+                std::cout << "Flow diff exceeded, flow_diff = " << rate_diff << std::endl;
-                return false;
+                ExceedInformation info;
+                info.group_name_ = name();
+                info.surplus_ = rate_diff;
+                info.well_index_ = self_index_;
+                result.fluid_rate_.push_back(info);
             }
             
         }
-        return true;
     }
 
     WellsGroupInterface* WellNode::findGroup(std::string name_of_node)
@@ -372,6 +374,9 @@ namespace Opm
             if (type == "GRUP") {
                 return ProductionSpecification::GRUP;
             }
+            if (type == "BHP") {
+                return ProductionSpecification::BHP;
+            }
 
             THROW("Unknown type " << type << ", could not convert to ControlMode.");
         }

opm/core/WellsGroup.hpp

@@ -86,7 +86,7 @@ namespace Opm
 
         void addChild(std::tr1::shared_ptr<WellsGroupInterface> child);
         
-        bool conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate,
+        void conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate,
                            const UnstructuredGrid& grid, const std::vector<double>& saturations, const struct Wells* wells,
                                    int index_of_well, WellControlResult& result, double epsilon = 1e-8);
         
@@ -108,7 +108,7 @@ namespace Opm
                 InjectionSpecification inj_spec);
 
         virtual WellsGroupInterface* findGroup(std::string name_of_node);
-        virtual bool conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
+        virtual void conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, 
                            const UnstructuredGrid& grid, const std::vector<double>& saturations, 
                            WellControlResult& result, double epsilon=1e-8);
         virtual bool isLeafNode() const;

opm/core/WellsManager.cpp

@@ -519,15 +519,7 @@ namespace Opm
         for (size_t i = 0; i < well_data.size(); i++) {
             if (well_collection_.getLeafNodes()[i]->prodSpec().control_mode_ == ProductionSpecification::GRUP) {
                 switch (well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_type_ ) {
-                case ProductionSpecification::OIL:
+                
-                {
-                    well_data[i].control = RATE;
-                    double parent_oil_rate = well_collection_.getLeafNodes()[i]->getParent()->prodSpec().oil_max_rate_;
-                    double guide_rate = well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_;
-                    well_data[i].target = guide_rate * parent_oil_rate;
-                    std::cout << "Applying guide rate" << std::endl;
-                    break;
-                }
                 case ProductionSpecification::NONE_GRT:
                 {
                     // Will use the group control type:
@@ -539,14 +531,20 @@ namespace Opm
                         well_data[i].target = guide_rate * parent_prod_spec.liquid_max_rate_;
                         well_data[i].control = RATE;
                         break;
+                    default:
+                        THROW("Unhandled group control mode " << parent_prod_spec.control_mode_);
+                        break;
                     }
                     
-                }   
+                }
+                default:
+                    break;
                 }
             }
 
             if (well_collection_.getLeafNodes()[i]->injSpec().control_mode_ == InjectionSpecification::GRUP) {
                 if (well_collection_.getLeafNodes()[i]->prodSpec().guide_rate_type_ == ProductionSpecification::RAT) {
+                    well_data[i].injected_phase = WATER; // Default for now.
                     well_data[i].control = RATE;
                     well_data[i].type = INJECTOR;
                     double parent_surface_rate = well_collection_.getLeafNodes()[i]->getParent()->injSpec().surface_flow_max_rate_;
@@ -558,7 +556,7 @@ namespace Opm
         }
 
 
-        
+        std::cout << "Making well structs" << std::endl;
 	// Set up the Wells struct.
 	w_ = wells_create(num_wells, num_perfs);
 	if (!w_) {
@@ -576,6 +574,11 @@ namespace Opm
 		wi[perf] = wellperf_data[w][perf].well_index;
 	    }
 	    const double* zfrac = (well_data[w].type == INJECTOR) ? fracs[well_data[w].injected_phase] : 0;
+            
+            // DIRTY DIRTY HACK
+            if(well_data[w].type == INJECTOR && well_data[w].injected_phase < 0 || well_data[w].injected_phase > 2){
+                zfrac = fracs[WATER];
+            }
 
             int ok = wells_add(well_data[w].type, well_data[w].reference_bhp_depth, nperf,
 			       zfrac, &cells[0], &wi[0], w_);
@@ -590,6 +593,7 @@ namespace Opm
 	    }
 	}
         
+        std::cout << "Made well struct" << std::endl;
         for(size_t i = 0; i < well_collection_.getLeafNodes().size(); i++) {
             WellNode* node = static_cast<WellNode*>(well_collection_.getLeafNodes()[i].get());
             
@@ -624,6 +628,120 @@ namespace Opm
         return well_collection_;
     }
 
+    
+    /// Apply control results
+    /// \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
+        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

@@ -20,6 +20,7 @@
 #ifndef OPM_WELLSMANAGER_HEADER_INCLUDED
 #define OPM_WELLSMANAGER_HEADER_INCLUDED
 #include <opm/core/WellCollection.hpp>
+#include <opm/core/WellsGroup.hpp>
 struct Wells;
 struct UnstructuredGrid;
 
@@ -49,6 +50,18 @@ namespace Opm
 
 	/// Destructor.
 	~WellsManager();
+        
+        /// Apply control results
+        /// \param[in] result The result of a run to conditionsMet on WellCollection
+        void applyControl(const WellControlResult& result);
+        
+        /// 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);
 
 	/// Access the managed Wells.
 	/// The method is named similarly to c_str() in std::string,