Merge pull request #627 from totto82/well_potentials

Compute well potentials

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -799,7 +799,7 @@ namespace detail {
 
     template <class Grid, class Implementation>
     void BlackoilModelBase<Grid, Implementation>::computeWellConnectionPressures(const SolutionState& state,
-                                                                        const WellState& xw)
+                                                                                 const WellState& xw)
     {
         if( ! localWellsActive() ) return ;
 

opm/autodiff/SimulatorBase.hpp

@@ -157,6 +157,12 @@ namespace Opm
                     const BlackoilState&            x,
                     WellState& xw);
 
+        void computeWellPotentials(const Wells*                    wells,
+                                   const BlackoilState&            x,
+                                   const WellState& xw,
+                                   std::vector<double>& well_potentials);
+
+
         // Data.
         typedef RateConverter::
         SurfaceToReservoirVoidage< BlackoilPropsAdInterface,

opm/autodiff/SimulatorBase_impl.hpp

@@ -1,6 +1,6 @@
 /*
   Copyright 2013 SINTEF ICT, Applied Mathematics.
-  Copyright 2014 IRIS AS
+  Copyright 2014-2016 IRIS AS
   Copyright 2015 Andreas Lauser
 
   This file is part of the Open Porous Media project (OPM).
@@ -22,6 +22,7 @@
 #include <algorithm>
 
 #include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
+#include <opm/core/well_controls.h>
 
 namespace Opm
 {
@@ -120,6 +121,8 @@ namespace Opm
         unsigned int totalNonlinearIterations = 0;
         unsigned int totalLinearIterations = 0;
 
+        std::vector<double> well_potentials;
+
         // Main simulation loop.
         while (!timer.done()) {
             // Report timestep.
@@ -139,7 +142,8 @@ namespace Opm
                                        Opm::UgGridHelpers::cell2Faces(grid_),
                                        Opm::UgGridHelpers::beginFaceCentroids(grid_),
                                        props_.permeability(),
-                                       is_parallel_run_);
+                                       is_parallel_run_,
+                                       well_potentials);
             const Wells* wells = wells_manager.c_wells();
             WellState well_state;
             well_state.init(wells, state, prev_well_state);
@@ -218,6 +222,10 @@ namespace Opm
             // Increment timer, remember well state.
             ++timer;
             prev_well_state = well_state;
+            // Compute Well potentials (only used to determine default guide rates for group controlled wells)
+            // TODO: add some logic to avoid unnecessary calulations of well potentials.
+            asImpl().computeWellPotentials(wells, state, well_state, well_potentials);
+
         }
 
         // Write final simulation state.
@@ -380,6 +388,43 @@ namespace Opm
         return std::unique_ptr<Solver>(new Solver(solver_param_, std::move(model)));
     }
 
+    template <class Implementation>
+    void SimulatorBase<Implementation>::computeWellPotentials(const Wells* wells,
+                                                              const BlackoilState& x,
+                                                              const WellState& xw,
+                                                              std::vector<double>& well_potentials)
+    {
+        const int nw = wells->number_of_wells;
+        const int np = wells->number_of_phases;
+
+        well_potentials.clear();
+        well_potentials.resize(nw*np,0.0);       
+        for (int w = 0; w < nw; ++w) {
+            for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                const double well_cell_pressure = x.pressure()[wells->well_cells[perf]];
+                const double drawdown_used = well_cell_pressure - xw.perfPress()[perf];
+                const WellControls* ctrl = wells->ctrls[w];
+                const int nwc = well_controls_get_num(ctrl);
+                //Loop over all controls until we find a BHP control
+                //that specifies what we need...
+                double bhp = 0.0;
+                for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
+                    if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
+                        bhp = well_controls_iget_target(ctrl, ctrl_index);
+                    }
+                    // TODO: do something for thp;
+                }
+                // Calculate the pressure difference in the well perforation
+                const double dp = xw.perfPress()[perf] - xw.bhp()[w];
+                const double drawdown_maximum = well_cell_pressure - (bhp + dp);
+
+                for (int phase = 0; phase < np; ++phase) {
+                    well_potentials[w*np + phase] += (drawdown_maximum / drawdown_used * xw.perfPhaseRates()[perf*np + phase]);
+                }
+            }
+        }
+    }
+
     template <class Implementation>
     void SimulatorBase<Implementation>::computeRESV(const std::size_t               step,
                                                     const Wells*                    wells,