merge.

opm/core/InjectionSpecification.cpp

@@ -3,11 +3,15 @@ namespace Opm
 {
 
     InjectionSpecification::InjectionSpecification()
-    : injector_type_(WATER), control_mode_(NONE), surface_flow_max_rate_(1e100),
-    reinjection_fraction_target_(0.0), BHP_limit_(1e100), fluid_volume_max_rate_(1e100)
+    : injector_type_(WATER),
+      control_mode_(NONE),
+      surface_flow_max_rate_(1e100),
+      reinjection_fraction_target_(0.0),
+      fluid_volume_max_rate_(1e100),
+      BHP_limit_(1e100)
     {
 
     }
 
 
-} // namespace Opm
+} // namespace Opm

opm/core/WellsGroup.cpp

@@ -7,6 +7,7 @@
 
 #include <opm/core/WellsGroup.hpp>
 #include <cmath>
+#include <opm/core/newwells.h>
 
 namespace Opm
 {
@@ -130,6 +131,7 @@ namespace Opm
         
         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]) {
@@ -145,6 +147,7 @@ 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;
         }
         }
@@ -153,12 +156,26 @@ namespace Opm
             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; 
+            
+            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) {
+                // 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;
         }
         return true;

opm/core/utility/miscUtilities.cpp

@@ -406,12 +406,12 @@ namespace Opm
     
     
     void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
-            const std::vector<double>& densities, std::vector<double>& wdp) 
+            const std::vector<double>& densities, std::vector<double>& wdp, bool per_grid_cell) 
     {
         const size_t np = densities.size();
-
+        const int nw = wells.number_of_wells;
         // Simple for now:
-        for(int i = 0; i < wells.number_of_wells; i++) {
+        for(int i = 0; i < nw; 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];
@@ -421,15 +421,25 @@ namespace Opm
 
                 double saturation_sum = 0.0;
                 for(size_t p = 0; p < np; p++) {
-                    saturation_sum += saturations[np*cell + p];
+                    if(per_grid_cell) {
+                        saturation_sum += saturations[i*nw*np + j*np + p];
+                    }
+                    else {
+                        saturation_sum += saturations[np*cell + p];
+                    }
                 }
                 if(saturation_sum == 0) {
                     saturation_sum = 1.0;
                 }
                 double density = 0.0;
                 for(size_t p = 0; p < np; p++) {
-                    // Is this a smart way of doing it?
-                    density += saturations[np*cell + p] * densities[p] / saturation_sum;
+                    if(per_grid_cell) {
+                        density += saturations[i*nw*np + j*np + p] * densities[p] / saturation_sum;
+                    }
+                    else {
+                        // Is this a smart way of doing it?
+                        density += saturations[np*cell + p] * densities[p] / saturation_sum;
+                    }
                 }
                 
                 // Is the sign correct?

opm/core/utility/miscUtilities.hpp

@@ -180,13 +180,15 @@ namespace Opm
     /// \param[in] wells        Wells that need their wdp calculated.
     /// \param[in] grid         The associated grid to make cell lookups.
     /// \param[in] saturations  A vector of weights for each cell for each phase 
-    ///                         in the grid. So for cell i, 
+    ///                         in the grid (or well, see per_grid_cell parameter). So for cell i, 
     ///                         saturations[i*densities.size() + p] should give the weight
     ///                         of phase p in cell i.
     /// \param[in] densities    Density for each phase.
     /// \param[out] wdp         Will contain, for each well, the wdp of the well.
+    /// \param[in] per_grid_cell Whether or not the saturations are per grid cell or per 
+    ///                          well cell.
     void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
-                    const std::vector<double>& densities, std::vector<double>& wdp);
+                    const std::vector<double>& densities, std::vector<double>& wdp, bool per_grid_cell = true);
     
     /// Computes (sums) the flow rate for each well. 
     /// \param[in] wells                The wells for which the flow rate should be computed.