Fixed bug in computeWDP(), add gravity argument. Make WellReport output in friendly units.

examples/spu_2p.cpp

@@ -142,6 +142,19 @@ static void outputWaterCut(const Opm::Watercut& watercut,
 }
 
 
+static void outputWellReport(const Opm::WellReport& wellreport,
+                             const std::string& output_dir)
+{
+    // Write well report.
+    std::string fname = output_dir  + "/wellreport.txt";
+    std::ofstream os(fname.c_str());
+    if (!os) {
+        THROW("Failed to open " << fname);
+    }
+    wellreport.write(os);
+}
+
+
 
 // --------------- Types needed to define transport solver ---------------
 
@@ -496,6 +509,15 @@ main(int argc, char** argv)
               << "    " << init_satvol[1]/tot_porevol_init << std::endl;
     Opm::Watercut watercut;
     watercut.push(0.0, 0.0, 0.0);
+    Opm::WellReport wellreport;
+    std::vector<double> well_bhp;
+    std::vector<double> well_perfrates;
+    if (wells->c_wells()) {
+        const int nw = wells->c_wells()->number_of_wells;
+        well_bhp.resize(nw, 0.0);
+        well_perfrates.resize(wells->c_wells()->well_connpos[nw], 0.0);
+        wellreport.push(*props, *wells->c_wells(), state.saturation(), 0.0, well_bhp, well_perfrates);
+    }
     for (; !simtimer.done(); ++simtimer) {
         // Report timestep and (optionally) write state to disk.
         simtimer.report(std::cout);
@@ -511,10 +533,8 @@ main(int argc, char** argv)
         }
         std::vector<double> wdp;
         if (wells->c_wells()) {
-            Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), wdp, true);
+            Opm::computeWDP(*wells->c_wells(), *grid->c_grid(), state.saturation(), props->density(), gravity[2], true, wdp);
         }
-        std::vector<double> well_bhp;
-        std::vector<double> well_perfrates;
         pressure_timer.start();
         if (rock_comp->isActive()) {
             rc.resize(num_cells);
@@ -639,6 +659,11 @@ main(int argc, char** argv)
         watercut.push(simtimer.currentTime() + simtimer.currentStepLength(),
                       produced[0]/(produced[0] + produced[1]),
                       tot_produced[0]/tot_porevol_init);
+        if (wells->c_wells()) {
+            wellreport.push(*props, *wells->c_wells(), state.saturation(),
+                            simtimer.currentTime() + simtimer.currentStepLength(),
+                            well_bhp, well_perfrates);
+        }
     }
     total_timer.stop();
 
@@ -650,6 +675,9 @@ main(int argc, char** argv)
     if (output) {
         outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir);
         outputWaterCut(watercut, output_dir);
+        if (wells->c_wells()) {
+            outputWellReport(wellreport, output_dir);
+        }
     }
 
     destroy_transport_source(tsrc);

examples/wells_example.cpp

@@ -57,7 +57,7 @@ int main(int argc, char** argv) {
     computeTotalMobilityOmega(incomp_properties, all_cells, state.saturation(), totmob, omega);
     
     std::vector<double> wdp;
-    computeWDP(*wells.c_wells(), *grid.c_grid(), state.saturation(), incomp_properties.density(), wdp, true);
+    computeWDP(*wells.c_wells(), *grid.c_grid(), state.saturation(), incomp_properties.density(), gravity[2], true, wdp);
     
     std::vector<double> src;
     Opm::FlowBCManager bcs;

opm/core/utility/miscUtilities.cpp

@@ -419,7 +419,8 @@ namespace Opm
 
 
     void computeWDP(const Wells& wells, const UnstructuredGrid& grid, const std::vector<double>& saturations,
-                    const double* densities, std::vector<double>& wdp, bool per_grid_cell)
+                    const double* densities, const double gravity, const bool per_grid_cell,
+                    std::vector<double>& wdp)
     {
         const int nw = wells.number_of_wells;
         const size_t np = per_grid_cell ?
@@ -456,7 +457,7 @@ namespace Opm
                 }
 
                 // Is the sign correct?
-                wdp.push_back(density * (cell_depth - depth_ref));
+                wdp.push_back(density * (cell_depth - depth_ref) * gravity);
             }
         }
     }
@@ -510,14 +511,14 @@ namespace Opm
         }
         const double* visc = props.viscosity();
         std::vector<double> data_now;
-        data_now.reserve(1 + 2*nw);
+        data_now.reserve(1 + 3*nw);
-        data_now.push_back(time);
+        data_now.push_back(time/unit::day);
         for (int w = 0; w < nw; ++w) {
-            data_now.push_back(well_bhp[w]);
+            data_now.push_back(well_bhp[w]/(unit::barsa));
             double well_rate_total = 0.0;
             double well_rate_water = 0.0;
             for (int perf = wells.well_connpos[w]; perf < wells.well_connpos[w + 1]; ++perf) {
-                const double perf_rate = well_perfrates[perf];
+                const double perf_rate = well_perfrates[perf]*(unit::day/unit::second);
                 well_rate_total += perf_rate;
                 if (perf_rate > 0.0) {
                     // Injection.
@@ -534,7 +535,11 @@ namespace Opm
                 }
             }
             data_now.push_back(well_rate_total);
-            data_now.push_back(well_rate_water/well_rate_total);
+            if (well_rate_total == 0.0) {
+                data_now.push_back(0.0);
+            } else {
+                data_now.push_back(well_rate_water/well_rate_total);
+            }
         }
         data_.push_back(data_now);
     }

opm/core/utility/miscUtilities.hpp

@@ -192,7 +192,8 @@ namespace Opm
     /// \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 double* densities, std::vector<double>& wdp, bool per_grid_cell);
+                    const double* densities, const double gravity, const bool per_grid_cell,
+                    std::vector<double>& wdp);
     
     /// Computes (sums) the flow rate for each well. 
     /// \param[in] wells                The wells for which the flow rate should be computed.