Added reading of WGRUPCON, also added guide_rates to the well groups. Fixed a small bug in calculateWDP

opm/core/WellCollection.cpp

@@ -99,4 +99,10 @@ namespace Opm
         
     }
 
+    void WellCollection::calculateGuideRates()
+    {
+        for(size_t i = 0; i < roots_.size(); i++) {
+            roots_[i]->calculateGuideRates();
+        }
+    }
 }

opm/core/WellCollection.hpp

@@ -43,6 +43,8 @@ namespace Opm
         bool conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate) const;
         
         const std::vector<std::tr1::shared_ptr<WellsGroupInterface> >& getLeafNodes() const;
+        
+        void calculateGuideRates();
     private:
         // To account for the possibility of a forest
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > roots_;

opm/core/WellsGroup.cpp

@@ -46,6 +46,32 @@ namespace Opm
         parent_ = parent;
     }
 
+    const ProductionSpecification& WellsGroupInterface::prodSpec() const
+    {
+        return production_specification_;
+    }
+
+    /// Injection specifications for the well or well group.
+
+    const InjectionSpecification& WellsGroupInterface::injSpec() const
+    {
+        return injection_specification_;
+    }
+
+    /// Production specifications for the well or well group.
+
+    ProductionSpecification& WellsGroupInterface::prodSpec()
+    {
+        return production_specification_;
+    }
+
+    /// Injection specifications for the well or well group.
+
+    InjectionSpecification& WellsGroupInterface::injSpec()
+    {
+        return injection_specification_;
+    }
+
     WellsGroupInterface* WellsGroup::findGroup(std::string name_of_node)
     {
         if (name() == name_of_node) {
@@ -64,6 +90,26 @@ namespace Opm
     }
 
     
+    void WellsGroup::calculateGuideRates()
+    {
+        double guide_rate_sum = 0.0;
+        for(size_t i = 0; i < children_.size(); i++) {
+            if(children_[i]->isLeafNode()) {
+                guide_rate_sum += children_[i]->prodSpec().guide_rate_;
+            }
+            else
+            {
+                children_[i]->calculateGuideRates();
+            }
+        }
+        if(guide_rate_sum != 0.0) {
+            for(size_t i = 0; i < children_.size(); i++) {
+                children_[i]->prodSpec().guide_rate_ /= guide_rate_sum;
+            }
+        }
+    }
+    
+    
     bool WellsGroup::conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate,
             const struct Wells* wells, int index_of_well)
     {
@@ -115,11 +161,17 @@ namespace Opm
         return true;
     }
 
-    void WellNode::setWellsPointer(const struct Wells* wells, int self_index) {
+    void WellNode::setWellsPointer(const struct Wells* wells, int self_index)
+    {
         wells_ = wells;
         self_index_ = self_index;
     }
     
+    void WellNode::calculateGuideRates()
+    {
+        // Empty
+    }
+
     namespace
     {
 
@@ -200,7 +252,6 @@ namespace Opm
             THROW("Unknown type " << type << ", could not convert to ControlMode.");
         }
 
-    
         ProductionSpecification::Procedure toProductionProcedure(std::string type)
         {
             if (type == "NONE") {

opm/core/WellsGroup.hpp

@@ -28,6 +28,12 @@ namespace Opm
         /// Injection specifications for the well or well group.
         const InjectionSpecification& injSpec() const;
         
+        /// Production specifications for the well or well group.
+        ProductionSpecification& prodSpec();
+        
+        /// Injection specifications for the well or well group.
+        InjectionSpecification& injSpec();
+        
         /// \returns true if the object is a leaf node (WellNode), false otherwise.
         virtual bool isLeafNode() const;
         
@@ -36,6 +42,8 @@ namespace Opm
         virtual WellsGroupInterface* findGroup(std::string name_of_node) = 0;
         
         void setParent(WellsGroupInterface* parent);
+        
+        virtual void calculateGuideRates() = 0;
     protected:
            WellsGroupInterface* parent_;
 
@@ -60,6 +68,8 @@ namespace Opm
         
         bool conditionsMet(const std::vector<double>& well_bhp, const std::vector<double>& well_rate, const struct Wells* wells,
                                    int index_of_well);
+        
+        virtual void calculateGuideRates();
     private:
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > children_;
     };
@@ -79,6 +89,7 @@ namespace Opm
         
         void setWellsPointer(const struct Wells* wells, int self_index);
         
+        virtual void calculateGuideRates();
     private:
         const struct Wells* wells_;
         int self_index_;

opm/core/WellsManager.cpp

@@ -208,6 +208,9 @@ namespace Opm
         std::vector<WellData> well_data;
 	std::vector<std::vector<PerfData> > wellperf_data;
         
+        // For easy lookup:
+        std::map<std::string, int> well_names_to_index;
+
 	// Get WELSPECS data
 	const WELSPECS& welspecs = deck.getWELSPECS();
 	const int num_wells = welspecs.welspecs.size();
@@ -218,6 +221,7 @@ namespace Opm
 	    well_names.push_back(welspecs.welspecs[w].name_);
 	    WellData wd;
 	    well_data.push_back(wd);
+            well_names_to_index[welspecs.welspecs[w].name_] = w;
             well_data.back().reference_bhp_depth = welspecs.welspecs[w].datum_depth_BHP_;
             if (welspecs.welspecs[w].datum_depth_BHP_ < 0.0) {
                 // Set refdepth to a marker value, will be changed
@@ -505,32 +509,44 @@ namespace Opm
 	    }
 	}
         
-        WellCollection wells_;
         if (deck.hasField("GRUPTREE")) {
             std::cout << "Found gruptree" << std::endl;
             const GRUPTREE& gruptree = deck.getGRUPTREE();
             
             std::map<std::string, std::string>::const_iterator it = gruptree.tree.begin();
             for( ; it != gruptree.tree.end(); ++it) {
-                wells_.addChild(it->first, it->second, deck);
+                well_collection_.addChild(it->first, it->second, deck);
             }
         }
         
         for (size_t i = 0; i < welspecs.welspecs.size(); ++i) {
             WelspecsLine line = welspecs.welspecs[i];
-            wells_.addChild(line.name_, line.group_, deck);
+            well_collection_.addChild(line.name_, line.group_, deck);
         }
         
-        for(size_t i = 0; i < wells_.getLeafNodes().size(); i++) {
+        for(size_t i = 0; i < well_collection_.getLeafNodes().size(); i++) {
-            WellNode* node = static_cast<WellNode*>(wells_.getLeafNodes()[i].get());
+            WellNode* node = static_cast<WellNode*>(well_collection_.getLeafNodes()[i].get());
             
             // We know that getLeafNodes() is ordered the same way as they're indexed in w_
             node->setWellsPointer(w_, i);
         }
         
         // Set the guide rates:
-        if(deck.hasField("WCONPROD")) {
+        if(deck.hasField("WGRUPCON")) {
-            WCONPROD wconprod = deck.getWCONPROD();
+            std::cout << "Found Wgrupcon" << std::endl;
+            WGRUPCON wgrupcon = deck.getWGRUPCON();
+            const std::vector<WgrupconLine>& lines = wgrupcon.wgrupcon;
+            std::cout << well_collection_.getLeafNodes().size() << std::endl;
+            for(size_t i = 0; i < lines.size(); i++) {
+                std::string name = lines[i].well_;
+                int index = well_names_to_index[name];
+                ASSERT(well_collection_.getLeafNodes()[index]->name() == name);
+                well_collection_.getLeafNodes()[index]->prodSpec().guide_rate_ = lines[i].guide_rate_;
+                well_collection_.getLeafNodes()[index]->prodSpec().guide_rate_type_ 
+                        = lines[i].phase_ == "OIL" ? ProductionSpecification::OIL : ProductionSpecification::RAT;
+            }
+            
+            well_collection_.calculateGuideRates();
         }
     }
 

opm/core/utility/miscUtilities.cpp

@@ -417,10 +417,17 @@ namespace Opm
                 // Is this correct wrt. depth_ref?
                 double cell_depth = grid.cell_centroids[3*cell+2];
 
+                double saturation_sum = 0.0;
+                for(size_t i = 0; i < densities.size(); i++) {
+                    saturation_sum += saturations[densities.size()*cell + i];
+                }
+                if(saturation_sum == 0) {
+                    saturation_sum = 1.0;
+                }
                 double density = 0.0;
                 for(size_t i = 0; i < densities.size(); i++) {
                     // Is this a smart way of doing it?
-                    density += saturations[densities.size()*cell+i]*densities[i];
+                    density += saturations[densities.size()*cell+i]*densities[i]/saturation_sum;
                 }
                 
                 // Is the sign correct?