Merge pull request #648 from totto82/fix_well_potentials

Fix well potential calculations

opm/autodiff/BlackoilModelBase.hpp

@@ -405,6 +405,13 @@ namespace Opm {
                     SolutionState& state,
                     WellState& well_state);
 
+        void
+        computeWellPotentials(const SolutionState& state,
+                              const std::vector<ADB>& mob_perfcells,
+                              const std::vector<ADB>& b_perfcells,
+                              WellState& well_state);
+
+
         void
         addWellFluxEq(const std::vector<ADB>& cq_s,
                       const SolutionState& state);

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -867,8 +867,9 @@ namespace detail {
         asImpl().addWellFluxEq(cq_s, state);
         asImpl().addWellContributionToMassBalanceEq(cq_s, state, well_state);
         asImpl().addWellControlEq(state, well_state, aliveWells);
-    }
+        asImpl().computeWellPotentials(state, mob_perfcells, b_perfcells, well_state);
 
+    }
 
 
 
@@ -1694,6 +1695,108 @@ namespace detail {
     }
 
 
+    template <class Grid, class Implementation>
+    void
+    BlackoilModelBase<Grid, Implementation>::
+    computeWellPotentials(const SolutionState& state,
+                          const std::vector<ADB>& mob_perfcells,
+                          const std::vector<ADB>& b_perfcells,
+                          WellState& well_state)
+    {
+        //only compute well potentials if they are needed
+        if (param_.compute_well_potentials_) {
+            const int nw = wells().number_of_wells;
+            const int np = wells().number_of_phases;
+            const Opm::PhaseUsage pu = fluid_.phaseUsage();
+            V bhps = V::Zero(nw);
+            for (int w = 0; w < nw; ++w) {
+                const WellControls* ctrl = wells().ctrls[w];
+                const int nwc = well_controls_get_num(ctrl);
+                //Loop over all controls until we find a BHP control
+                //or a THP control that specifies what we need.
+                //Pick the value that gives the most restrictive flow
+                for (int ctrl_index=0; ctrl_index < nwc; ++ctrl_index) {
+
+                    if (well_controls_iget_type(ctrl, ctrl_index) == BHP) {
+                        bhps[w] = well_controls_iget_target(ctrl, ctrl_index);
+                    }
+
+                    if (well_controls_iget_type(ctrl, ctrl_index) == THP) {
+                        double aqua = 0.0;
+                        double liquid = 0.0;
+                        double vapour = 0.0;
+
+                        if (active_[ Water ]) {
+                            aqua = well_state.wellRates()[w*np + pu.phase_pos[ Water ] ];
+                        }
+                        if (active_[ Oil ]) {
+                            liquid = well_state.wellRates()[w*np + pu.phase_pos[ Oil ] ];
+                        }
+                        if (active_[ Gas ]) {
+                            vapour = well_state.wellRates()[w*np + pu.phase_pos[ Gas ] ];
+                        }
+
+                        const int vfp        = well_controls_iget_vfp(ctrl, ctrl_index);
+                        const double& thp    = well_controls_iget_target(ctrl, ctrl_index);
+                        const double& alq    = well_controls_iget_alq(ctrl, ctrl_index);
+
+                        //Set *BHP* target by calculating bhp from THP
+                        const WellType& well_type = wells().type[w];
+
+                        const double gravity = detail::getGravity(geo_.gravity(), UgGridHelpers::dimensions(grid_));
+
+                        if (well_type == INJECTOR) {
+                            double dp = wellhelpers::computeHydrostaticCorrection(
+                                        wells(), w, vfp_properties_.getInj()->getTable(vfp)->getDatumDepth(),
+                                        stdWells().wellPerforationDensities(), gravity);
+                            const double bhp = vfp_properties_.getInj()->bhp(vfp, aqua, liquid, vapour, thp) - dp;
+                            // apply the strictest of the bhp controlls i.e. smallest bhp for injectors
+                            if ( bhp < bhps[w]) {
+                                bhps[w] = bhp;
+                            }
+                        }
+                        else if (well_type == PRODUCER) {
+                            double dp = wellhelpers::computeHydrostaticCorrection(
+                                        wells(), w, vfp_properties_.getProd()->getTable(vfp)->getDatumDepth(),
+                                        stdWells().wellPerforationDensities(), gravity);
+
+                            const double bhp = vfp_properties_.getProd()->bhp(vfp, aqua, liquid, vapour, thp, alq) - dp;
+                            // apply the strictest of the bhp controlls i.e. largest bhp for producers
+                            if ( bhp > bhps[w]) {
+                                bhps[w] = bhp;
+                            }
+                        }
+                        else {
+                            OPM_THROW(std::logic_error, "Expected PRODUCER or INJECTOR type of well");
+                        }
+                    }
+                }
+
+            }
+
+            // use bhp limit from control
+            SolutionState state0 = state;
+            asImpl().makeConstantState(state0);
+            state0.bhp = ADB::constant(bhps);
+
+            // compute well potentials
+            V aliveWells;
+            std::vector<ADB> well_potentials;
+            asImpl().stdWells().computeWellFlux(state0, fluid_.phaseUsage(), active_, mob_perfcells,  b_perfcells, aliveWells, well_potentials);
+
+            // store well potentials in the well state
+            // transform to a single vector instead of separate vectors pr phase
+            const int nperf = wells().well_connpos[nw];
+            V cq = superset(well_potentials[0].value(), Span(nperf, np, 0), nperf*np);
+            for (int phase = 1; phase < np; ++phase) {
+                cq += superset(well_potentials[phase].value(), Span(nperf, np, phase), nperf*np);
+            }
+            well_state.wellPotentials().assign(cq.data(), cq.data() + nperf*np);
+        }
+
+    }
+
+
 
 
 

opm/autodiff/BlackoilModelParameters.cpp

@@ -48,6 +48,7 @@ namespace Opm
         tolerance_wells_ = param.getDefault("tolerance_wells", tolerance_wells_ );
         solve_welleq_initially_ = param.getDefault("solve_welleq_initially",solve_welleq_initially_);
         update_equations_scaling_ = param.getDefault("update_equations_scaling", update_equations_scaling_);
+        compute_well_potentials_ = param.getDefault("compute_well_potentials", compute_well_potentials_);
     }
 
 
@@ -65,6 +66,7 @@ namespace Opm
         tolerance_wells_ = 1.0e-3;
         solve_welleq_initially_ = true;
         update_equations_scaling_ = false;
+        compute_well_potentials_ = false;
     }
 
 

opm/autodiff/BlackoilModelParameters.hpp

@@ -49,6 +49,10 @@ namespace Opm
         /// Update scaling factors for mass balance equations
         bool update_equations_scaling_;
 
+        /// Compute well potentials, needed to calculate default guide rates for group
+        /// controlled wells
+        bool compute_well_potentials_;
+
         /// Construct from user parameters or defaults.
         explicit BlackoilModelParameters( const parameter::ParameterGroup& param );
 

opm/autodiff/SimulatorBase_impl.hpp

@@ -120,7 +120,7 @@ namespace Opm
 
         unsigned int totalNonlinearIterations = 0;
         unsigned int totalLinearIterations = 0;
-
+        bool is_well_potentials_computed = param_.getDefault("compute_well_potentials", false );
         std::vector<double> well_potentials;
 
         // Main simulation loop.
@@ -222,12 +222,13 @@ 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);
+            // The well potentials are only computed if they are needed
+            // For now thay are only used to determine default guide rates for group controlled wells
+            if ( is_well_potentials_computed ) {
+                asImpl().computeWellPotentials(wells, state, well_state, well_potentials);
+            }
 
         }
-
         // Write final simulation state.
         output_writer_.writeTimeStep( timer, state, prev_well_state );
 
@@ -396,30 +397,12 @@ namespace Opm
     {
         const int nw = wells->number_of_wells;
         const int np = wells->number_of_phases;
-
         well_potentials.clear();
-        well_potentials.resize(nw*np,0.0);       
+        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]);
+                    well_potentials[w*np + phase] += xw.wellPotentials()[perf*np + phase];
                 }
             }
         }

opm/autodiff/WellStateFullyImplicitBlackoil.hpp

@@ -102,6 +102,9 @@ namespace Opm
                 current_controls_[w] = well_controls_get_current(wells->ctrls[w]);
             }
 
+            well_potentials_.clear();
+            well_potentials_.resize(nperf * np, 0.0);
+
             // intialize wells that have been there before
             // order may change so the mapping is based on the well name
             if( ! prevState.wellMap().empty() )
@@ -184,9 +187,14 @@ namespace Opm
         std::vector<int>& currentControls() { return current_controls_; }
         const std::vector<int>& currentControls() const { return current_controls_; }
 
+        /// One rate per phase and well connection.
+        std::vector<double>& wellPotentials() { return well_potentials_; }
+        const std::vector<double>& wellPotentials() const { return well_potentials_; }
+
     private:
         std::vector<double> perfphaserates_;
         std::vector<int> current_controls_;
+        std::vector<double> well_potentials_;
     };
 
 } // namespace Opm