Made a computeWDP-function

examples/wells_example.cpp

@@ -70,7 +70,8 @@ int main(int argc, char** argv) {
 
     std::vector<double> well_bhp;
     std::vector<double> well_rate;
-    //pressure_solver.solve(totmob, omega, src, bcs.c_bcs(), pressure, face_flux, well_bhp, well_rate);
+    pressure_solver.solve(totmob, omega, src, bcs.c_bcs(), pressure, face_flux, well_bhp, well_rate);
+    std::cout << "Solved" << std::endl;
     if(wells.wellCollection().conditionsMet(well_bhp, well_rate)) {
         std::cout << "Conditions met for wells!" << std::endl;
     }

opm/core/pressure/IncompTpfa.cpp

@@ -124,6 +124,7 @@ namespace Opm
         ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
         if (! src.empty()) { F.src = &src[0]; }
         F.bc = bcs;
+        F.totmob = &totmob[0];
 
 	ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
 
@@ -199,7 +200,7 @@ namespace Opm
         ifs_tpfa_forces F = { NULL, NULL, wells_, NULL, NULL };
         if (! src.empty()) { F.src = &src[0]; }
         F.bc = bcs;
-
+        F.totmob = &totmob[0];
 	ifs_tpfa_assemble(gg, &F, &trans_[0], &gpress_omegaweighted_[0], h_);
 
         // TODO: this is a hack, it would be better to handle this in a

opm/core/utility/miscUtilities.cpp

@@ -403,6 +403,31 @@ namespace Opm
                << data_[3*i+2] << '\n';
         }
     }
+    
+    
+    void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
+            const std::vector<double> densities, std::vector<double>& wdp) 
+    {
+        // Simple for now:
+        for(int i = 0; i < wells.number_of_wells; i++) {
+            double depth_ref = wells.depth_ref[i];
+            for(int j = wells.well_connpos[i]; j < wells.well_connpos[i+1]; j++) {
+                int cell = wells.well_cells[j];
+                
+                // Is this correct wrt. depth_ref?
+                double cell_depth = grid.cell_centroids[3*cell+2];
+                
+                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];
+                }
+                
+                // Is the sign correct?
+                wdp.push_back(density*(cell_depth-depth_ref));
+            }
+        }
+    }
 
 
 } // namespace Opm

opm/core/utility/miscUtilities.hpp

@@ -175,7 +175,11 @@ namespace Opm
     /// For this to be valid, the wells must be all rate-controlled and
     /// single-perforation.
     void wellsToSrc(const Wells& wells, const int num_cells, std::vector<double>& src);
+    
 
+    void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
+                    const std::vector<double>& densities, std::vector<double>& wdp);
+    
     /// Encapsulates the watercut curves.
     class Watercut
     {