Merged.

examples/wells_example.cpp

@@ -3,8 +3,11 @@
 #include <opm/core/utility/parameters/ParameterGroup.hpp>
 #include <opm/core/WellsManager.hpp>
 #include <opm/core/GridManager.hpp>
-
+#include <opm/core/pressure/IncompTpfa.hpp>
+#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
+#include <opm/core/linalg/LinearSolverUmfpack.hpp>
 #include "opm/core/newwells.h"
+#include "opm/core/grid.h"
 
 int main(int argc, char** argv) {
 
@@ -12,7 +15,7 @@ int main(int argc, char** argv) {
     using namespace Opm;
     ParameterGroup parameters( argc, argv, false );
     std::string file_name = parameters.getDefault<std::string>("inputdeck", "data.data");
-    
+
     // Read input file
     EclipseGridParser parser(file_name);
     std::cout << "Done!" << std::endl;
@@ -22,6 +25,18 @@ int main(int argc, char** argv) {
     // Finally handle the wells
     WellsManager wells(parser, *grid.c_grid(), NULL);
     
+    std::vector<int> global_cells(grid.c_grid()->global_cell, grid.c_grid()->global_cell + grid.c_grid()->number_of_cells);
+
+    std::cout << "ahoi" << std::endl;
+    double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
+    IncompPropertiesFromDeck incomp_properties(parser, global_cells);
+        std::cout << "there" << std::endl;
+
+    LinearSolverUmfpack umfpack_solver;
+    std::cout << "here" << std::endl;
+    IncompTpfa pressure_solver(*grid.c_grid(), incomp_properties.permeability(), 
+                               gravity, umfpack_solver,  wells.c_wells());
+    
     return 0;
 }
 

opm/core/WellCollection.cpp

@@ -41,11 +41,11 @@ namespace Opm
         }
         std::tr1::shared_ptr<WellsGroupInterface> child;
 
-        for (int i = 0; i < roots_.size(); ++i) {
+        for (size_t i = 0; i < roots_.size(); ++i) {
             if (roots_[i]->name() == child_name) {
                 child = roots_[i];
                 // We've found a new parent to the previously thought root, need to remove it
-                for(int j = i; j < roots_.size() - 1; ++j) {
+                for(size_t j = i; j < roots_.size() - 1; ++j) {
                     roots_[j] = roots_[j+1];
                 }
                 
@@ -63,13 +63,19 @@ namespace Opm
         }
         parent_as_group->addChild(child);
 
-
+        if(child->isLeafNode()) {
+            leaf_nodes_.push_back(child);
+        }
+    }
+    
+    const std::vector<std::tr1::shared_ptr<WellsGroupInterface> >& WellCollection::getLeafNodes() const {
+        return leaf_nodes_;
     }
 
     WellsGroupInterface* WellCollection::findNode(std::string name)
     {
 
-        for (int i = 0; i < roots_.size(); i++) {
+        for (size_t i = 0; i < roots_.size(); i++) {
             WellsGroupInterface* result = roots_[i]->findGroup(name);
             if (result) {
                 return result;

opm/core/WellCollection.hpp

@@ -38,9 +38,14 @@ namespace Opm
         void addChild(std::string child, std::string parent,
                 const EclipseGridParser& deck);
         
+        const std::vector<std::tr1::shared_ptr<WellsGroupInterface> >& getLeafNodes() const;
     private:
         // To account for the possibility of a forest
         std::vector<std::tr1::shared_ptr<WellsGroupInterface> > roots_;
+        
+        // This will be used to traverse the bottom nodes.
+        std::vector<std::tr1::shared_ptr<WellsGroupInterface> > leaf_nodes_;
+        
         WellsGroupInterface* findNode(std::string name);
     };
 

opm/core/WellsGroup.cpp

@@ -31,13 +31,17 @@ namespace Opm
     : WellsGroupInterface(name, prod_spec, inj_spec)
     {
     }
+    
+    bool WellsGroupInterface::isLeafNode() const {
+        return false;
+    }
 
     WellsGroupInterface* WellsGroup::findGroup(std::string name_of_node)
     {
         if (name() == name_of_node) {
             return this;
         } else {
-            for (int i = 0; i < children_.size(); i++) {
+            for (size_t i = 0; i < children_.size(); i++) {
                 WellsGroupInterface* result = children_[i]->findGroup(name_of_node);
                 if (result) {
                     return result;
@@ -68,6 +72,10 @@ namespace Opm
         }
 
     }
+    
+    bool WellNode::isLeafNode() const {
+        return true;
+    }
 
     surface_component toSurfaceComponent(std::string type)
     {
@@ -171,7 +179,7 @@ namespace Opm
         bool isWell = false;
         if (deck.hasField("WELSPECS")) {
             WELSPECS wspecs = deck.getWELSPECS();
-            for (int i = 0; i < wspecs.welspecs.size(); i++) {
+            for (size_t i = 0; i < wspecs.welspecs.size(); i++) {
                 if (wspecs.welspecs[i].name_ == name) {
                     isWell = true;
                     break;
@@ -185,7 +193,7 @@ namespace Opm
             InjectionSpecification injection_specification;
             if (deck.hasField("WCONINJE")) {
                 WCONINJE wconinje = deck.getWCONINJE();
-                for (int i = 0; i < wconinje.wconinje.size(); i++) {
+                for (size_t i = 0; i < wconinje.wconinje.size(); i++) {
                     if (wconinje.wconinje[i].well_ == name) {
                         WconinjeLine line = wconinje.wconinje[i];
                         injection_specification.BHP_limit_ = line.BHP_limit_;
@@ -199,7 +207,7 @@ namespace Opm
             ProductionSpecification production_specification;
             if (deck.hasField("WCONPROD")) {
                 WCONPROD wconprod = deck.getWCONPROD();
-                for (int i = 0; i < wconprod.wconprod.size(); i++) {
+                for (size_t i = 0; i < wconprod.wconprod.size(); i++) {
                     if (wconprod.wconprod[i].well_ == name) {
                         WconprodLine line = wconprod.wconprod[i];
                         production_specification.BHP_limit_ = line.BHP_limit_;
@@ -216,7 +224,7 @@ namespace Opm
             InjectionSpecification injection_specification;
             if (deck.hasField("GCONINJE")) {
                 GCONINJE gconinje = deck.getGCONINJE();
-                for (int i = 0; i < gconinje.gconinje.size(); i++) {
+                for (size_t i = 0; i < gconinje.gconinje.size(); i++) {
                     if (gconinje.gconinje[i].group_ == name) {
                         GconinjeLine line = gconinje.gconinje[i];
                         injection_specification.injector_type_ = toSurfaceComponent(line.injector_type_);
@@ -229,7 +237,7 @@ namespace Opm
             ProductionSpecification production_specification;
             if (deck.hasField("GCONPROD")) {
                 GCONPROD gconprod = deck.getGCONPROD();
-                for (int i = 0; i < gconprod.gconprod.size(); i++) {
+                for (size_t i = 0; i < gconprod.gconprod.size(); i++) {
                     if (gconprod.gconprod[i].group_ == name) {
                         GconprodLine line = gconprod.gconprod[i];
                         production_specification.oil_max_rate_ = line.oil_max_rate_;

opm/core/WellsGroup.hpp

@@ -17,9 +17,17 @@ namespace Opm
                 InjectionSpecification inj_spec);
         virtual ~WellsGroupInterface();
 
+        /// The unique identifier for the well or well group.
         const std::string& name();
+        
+        /// Production specifications for the well or well group.
         const ProductionSpecification& prodSpec() const;
+        
+        /// Injection specifications for the well or well group.
         const InjectionSpecification& injSpec() const;
+        
+        /// \returns true if the object is a leaf node (WellNode), false otherwise.
+        virtual bool isLeafNode() const;
 
         /// \returns the pointer to the WellsGroupInterface with the given name. NULL if 
         ///          the name is not found.a
@@ -51,6 +59,8 @@ namespace Opm
                 InjectionSpecification inj_spec);
 
         virtual WellsGroupInterface* findGroup(std::string name_of_node);
+        
+        virtual bool isLeafNode() const;
 
     };
 

opm/core/WellsManager.cpp

@@ -505,22 +505,22 @@ namespace Opm
 	    }
 	}
         
-        WellCollection well_collection;
+        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) {
-                well_collection.addChild(it->first, it->second, deck);
+                wells_.addChild(it->first, it->second, deck);
             }
         }
         
         if(deck.hasField("WELSPECS")) {
             WELSPECS welspecs = deck.getWELSPECS();
-            for(int i = 0; i < welspecs.welspecs.size(); ++i) {
+            for(size_t i = 0; i < welspecs.welspecs.size(); ++i) {
                 WelspecsLine line = welspecs.welspecs[i];
-                well_collection.addChild(line.name_, line.group_, deck);
+                wells_.addChild(line.name_, line.group_, deck);
             }
         }
         

opm/core/WellsManager.hpp

@@ -19,7 +19,7 @@
 
 #ifndef OPM_WELLSMANAGER_HEADER_INCLUDED
 #define OPM_WELLSMANAGER_HEADER_INCLUDED
-
+#include <opm/core/WellCollection.hpp>
 struct Wells;
 struct UnstructuredGrid;
 
@@ -61,6 +61,8 @@ namespace Opm
 	WellsManager& operator=(const WellsManager& other);
 	// The managed Wells.
 	Wells* w_;
+        
+        WellCollection well_collection_;
 
     };
 

opm/core/pressure/IncompTpfa.cpp

@@ -38,14 +38,18 @@ namespace Opm
     ///                          and N == g.number_of_cells.
     /// \param[in] gravity       Gravity vector. If nonzero, the array should
     ///                          have D elements.
+    /// \param[in] wells         The wells argument. Will be used in solution, 
+    ///                          is ignored if NULL
     IncompTpfa::IncompTpfa(const UnstructuredGrid& g,
 			   const double* permeability,
 			   const double* gravity,
-                           const LinearSolverInterface& linsolver)
+                           const LinearSolverInterface& linsolver,
+                           const struct Wells* wells)
 	: grid_(g),
           linsolver_(linsolver),
 	  htrans_(g.cell_facepos[ g.number_of_cells ]),
-	  trans_ (g.number_of_faces)
+	  trans_ (g.number_of_faces),
+          wells_(wells)
     {
 	UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
 	tpfa_htrans_compute(gg, permeability, &htrans_[0]);
@@ -110,7 +114,7 @@ namespace Opm
 	    }
 	}
 
-        ifs_tpfa_forces F = { NULL, NULL, NULL, NULL, NULL };
+        ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
         if (! src.empty()) { F.src = &src[0]; }
         F.bc = bcs;
 
@@ -172,7 +176,7 @@ namespace Opm
 	    }
 	}
 
-        ifs_tpfa_forces F = { NULL, NULL, NULL, NULL, NULL };
+        ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
         if (! src.empty()) { F.src = &src[0]; }
         F.bc = bcs;
 

opm/core/pressure/IncompTpfa.hpp

@@ -25,6 +25,7 @@
 
 struct UnstructuredGrid;
 struct ifs_tpfa_data;
+struct Wells;
 struct FlowBoundaryConditions;
 
 namespace Opm
@@ -47,10 +48,13 @@ namespace Opm
 	///                          and N == g.number_of_cells.
 	/// \param[in] gravity       Gravity vector. If nonzero, the array should
 	///                          have D elements.
+        /// \param[in] wells         The wells argument. Will be used in solution, 
+        ///                          is ignored if NULL
 	IncompTpfa(const UnstructuredGrid& g,
 		   const double* permeability,
 		   const double* gravity,
-                   const LinearSolverInterface& linsolver);
+                   const LinearSolverInterface& linsolver,
+                   const struct Wells* wells = 0);
 
 	/// Destructor.
 	~IncompTpfa();
@@ -113,7 +117,8 @@ namespace Opm
 	::std::vector<double> trans_ ;
 	::std::vector<double> gpress_;
 	::std::vector<double> gpress_omegaweighted_;
-
+        
+        const struct Wells* wells_;
 	struct ifs_tpfa_data* h_;
     };
 

opm/core/utility/writeVtkData.cpp

@@ -22,6 +22,7 @@
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <boost/lexical_cast.hpp>
 #include <set>
+#include <cmath>
 #include <algorithm>
 #include <iostream>
 #include <iterator>