Merge remote-tracking branch 'upstream/master'

opm/core/pressure/CompressibleTpfa.cpp

@@ -429,7 +429,7 @@ namespace Opm
             for (int phase = 0; phase < np; ++phase) {
                 int upwindc = -1;
                 if (c[0] >=0 && c[1] >= 0) {
-                    upwindc = (pot[0] < pot[1]) ? c[1] : c[0];
+                    upwindc = (pot[0][phase] < pot[1][phase]) ? c[1] : c[0];
                 } else {
                     upwindc = (c[0] >= 0) ? c[0] : c[1];
                 }

opm/core/props/rock/RockCompressibility.cpp

@@ -42,12 +42,16 @@ namespace Opm
     {
         if (deck.hasField("ROCKTAB")) {
             const table_t& rt = deck.getROCKTAB().rocktab_;
-            int n = rt[0][0].size();
+            if (rt.size() != 1) {
+                THROW("Can only handle a single region in ROCKTAB.");
+            }
+            const int n = rt[0][0].size();
             p_.resize(n);
             poromult_.resize(n);
             for (int i = 0; i < n; ++i) {
                 p_[i] = rt[0][0][i];
                 poromult_[i] = rt[0][1][i];
+                transmult_[i] = rt[0][2][i];
             }
         } else if (deck.hasField("ROCK")) {
             const ROCK& r = deck.getROCK();
@@ -70,11 +74,39 @@ namespace Opm
             const double cpnorm = rock_comp_*(pressure - pref_);
             return (1.0 + cpnorm + 0.5*cpnorm*cpnorm);
         } else {
-            // return Opm::linearInterpolation(p_, poromult_, pressure);
             return Opm::linearInterpolation(p_, poromult_, pressure);
         }
     }
 
+    double RockCompressibility::poroMultDeriv(double pressure) const
+    {
+        if (p_.empty()) {
+            // Approximating poro multiplier with a quadratic curve,
+            // we must use its derivative.
+            return rock_comp_ + 2 * rock_comp_ * rock_comp_ * (pressure - pref_);
+        } else {
+            return Opm::linearInterpolationDerivative(p_, poromult_, pressure);
+        }
+    }
+
+    double RockCompressibility::transMult(double pressure) const
+    {
+        if (p_.empty()) {
+            return 1.0;
+        } else {
+            return Opm::linearInterpolation(p_, transmult_, pressure);
+        }
+    }
+
+    double RockCompressibility::transMultDeriv(double pressure) const
+    {
+        if (p_.empty()) {
+            return 0.0;
+        } else {
+            return Opm::linearInterpolationDerivative(p_, transmult_, pressure);
+        }
+    }
+
     double RockCompressibility::rockComp(double pressure) const
     {
         if (p_.empty()) {

opm/core/props/rock/RockCompressibility.hpp

@@ -47,12 +47,22 @@ namespace Opm
         /// Porosity multiplier.
         double poroMult(double pressure) const;
 
+        /// Derivative of porosity multiplier with respect to pressure.
+        double poroMultDeriv(double pressure) const;
+
+        /// Transmissibility multiplier.
+        double transMult(double pressure) const;
+
+        /// Derivative of transmissibility multiplier with respect to pressure.
+        double transMultDeriv(double pressure) const;
+
         /// Rock compressibility = (d poro / d p)*(1 / poro).
         double rockComp(double pressure) const;
 
     private:
         std::vector<double> p_;
         std::vector<double> poromult_;
+        std::vector<double> transmult_;
         double pref_;
         double rock_comp_;
     };

opm/core/simulator/BlackoilState.hpp

@@ -22,6 +22,7 @@
 
 #include <opm/core/grid.h>
 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/utility/ErrorMacros.hpp>
 #include <vector>
 
 namespace Opm
@@ -47,6 +48,7 @@ namespace Opm
                 // but use available phase information instead.
                 sat_[num_phases*cell + 1] = 1.0;
             }
+            gor_.resize(g.number_of_cells, 0.0);
         }
 
         enum ExtremalSat { MinSat, MaxSat };
@@ -59,7 +61,11 @@ namespace Opm
                          const Opm::BlackoilPropertiesInterface& props,
                          ExtremalSat es)
         {
+            if (cells.empty()) {
+                return;
+            }
             const int n = cells.size();
+            ASSERT(n > 0);
             std::vector<double> smin(num_phases_*n);
             std::vector<double> smax(num_phases_*n);
             props.satRange(n, &cells[0], &smin[0], &smax[0]);
@@ -81,12 +87,14 @@ namespace Opm
         std::vector<double>& faceflux    () { return flux_  ; }
         std::vector<double>& surfacevol  () { return surfvol_; }
         std::vector<double>& saturation  () { return sat_   ; }
+        std::vector<double>& gasoilratio () { return gor_   ; }
 
         const std::vector<double>& pressure    () const { return press_ ; }
         const std::vector<double>& facepressure() const { return fpress_; }
         const std::vector<double>& faceflux    () const { return flux_  ; }
         const std::vector<double>& surfacevol  () const { return surfvol_; }
         const std::vector<double>& saturation  () const { return sat_   ; }
+        const std::vector<double>& gasoilratio () const { return gor_   ; }
 
     private:
         int num_phases_;
@@ -95,6 +103,7 @@ namespace Opm
         std::vector<double> flux_  ;
         std::vector<double> surfvol_;
         std::vector<double> sat_   ;
+        std::vector<double> gor_   ;
     };
 
 } // namespace Opm

opm/core/simulator/WellState.hpp

@@ -31,27 +31,46 @@ namespace Opm
     {
     public:
         /// Allocate and initialize if wells is non-null.
+        /// Also tries to give useful initial values to the bhp() and
+        /// wellRates() fields, depending on controls.  The
+        /// perfRates() field is filled with zero, and perfPress()
+        /// with -1e100.
         template <class State>
         void init(const Wells* wells, const State& state)
         {
             if (wells) {
                 const int nw = wells->number_of_wells;
+                const int np = wells->number_of_phases;
                 bhp_.resize(nw);
-                // Initialize bhp to be target pressure
+                wellrates_.resize(nw * np, 0.0);
-                // if bhp-controlled well, otherwise set
-                // to pressure in first perforation cell.
                 for (int w = 0; w < nw; ++w) {
                     const WellControls* ctrl = wells->ctrls[w];
-
+                    // Initialize bhp to be target pressure if
+                    // bhp-controlled well, otherwise set to a little
+                    // above or below (depending on if the well is an
+                    // injector or producer) pressure in first perforation
+                    // cell.
                     if ((ctrl->current < 0) || // SHUT
                         (ctrl->type[ctrl->current] != BHP)) {
-                        const int cell = wells->well_cells[wells->well_connpos[w]];
+                        const int first_cell = wells->well_cells[wells->well_connpos[w]];
-                        bhp_[w] = state.pressure()[cell];
+                        const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
-                    }
+                        bhp_[w] = safety_factor*state.pressure()[first_cell];
-                    else {
+                    } else {
                         bhp_[w] = ctrl->target[ctrl->current];
                     }
+                    // Initialize well rates to match controls if type is SURFACE_RATE
+                    if ((ctrl->current >= 0) && // open well
+                        (ctrl->type[ctrl->current] == SURFACE_RATE)) {
+                        const double rate_target = ctrl->target[ctrl->current];
+                        for (int p = 0; p < np; ++p) {
+                            const double phase_distr = ctrl->distr[np * ctrl->current + p];
+                            wellrates_[np*w + p] = rate_target * phase_distr;
                         }
+                    }
+                }
+                // The perforation rates and perforation pressures are
+                // not expected to be consistent with bhp_ and wellrates_
+                // after init().
                 perfrates_.resize(wells->well_connpos[nw], 0.0);
                 perfpress_.resize(wells->well_connpos[nw], -1e100);
             }
@@ -61,6 +80,10 @@ namespace Opm
         std::vector<double>& bhp() { return bhp_; }
         const std::vector<double>& bhp() const { return bhp_; }
 
+        /// One rate per well and phase.
+        std::vector<double>& wellRates() { return wellrates_; }
+        const std::vector<double>& wellRates() const { return wellrates_; }
+
         /// One rate per well connection.
         std::vector<double>& perfRates() { return perfrates_; }
         const std::vector<double>& perfRates() const { return perfrates_; }
@@ -71,6 +94,7 @@ namespace Opm
 
     private:
         std::vector<double> bhp_;
+        std::vector<double> wellrates_;
         std::vector<double> perfrates_;
         std::vector<double> perfpress_;
     };

opm/core/simulator/initState_impl.hpp

@@ -599,6 +599,29 @@ namespace Opm
     }
 
 
+    /// Initialize a blackoil state from input deck.
+    template <class Props, class State>
+    void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
+                                   const Props& props,
+                                   const EclipseGridParser& deck,
+                                   const double gravity,
+                                   State& state)
+    {
+        initStateFromDeck(grid, props, deck, gravity, state);
+        initBlackoilSurfvol(grid, props, state);
+        if (deck.hasField("RS")) {
+            const std::vector<double>& rs_deck = deck.getFloatingPointValue("RS");
+            const int num_cells = grid.number_of_cells;
+            for (int c = 0; c < num_cells; ++c) {
+                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                state.gasoilratio()[c] = rs_deck[c_deck];
+            }
+        } else {
+            THROW("Temporarily, we require the RS field.");
+        }
+    }
+
+
 } // namespace Opm
 
 

opm/core/wells/WellsManager.cpp

@@ -505,13 +505,33 @@ namespace Opm
                         int control_pos[5] = { -1, -1, -1, -1, -1 };
                         if (ok && wci_line.surface_flow_max_rate_ >= 0.0) {
                             control_pos[InjectionControl::RATE] = w_->ctrls[wix]->num;
-                            const double distr[3] = { 1.0, 1.0, 1.0 };
+                            double distr[3] = { 0.0, 0.0, 0.0 };
+                            if (wci_line.injector_type_ == "WATER") {
+                                distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                            } else if (wci_line.injector_type_ == "OIL") {
+                                distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                            } else if (wci_line.injector_type_ == "GAS") {
+                                distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                            } else {
+                                THROW("Injector type " << wci_line.injector_type_ << " not supported."
+                                      "WellsManager only supports WATER, OIL and GAS injector types.");
+                            }
                             ok = append_well_controls(SURFACE_RATE, wci_line.surface_flow_max_rate_,
                                                       distr, wix, w_);
                         }
                         if (ok && wci_line.reservoir_flow_max_rate_ >= 0.0) {
                             control_pos[InjectionControl::RESV] = w_->ctrls[wix]->num;
-                            const double distr[3] = { 1.0, 1.0, 1.0 };
+                            double distr[3] = { 0.0, 0.0, 0.0 };
+                            if (wci_line.injector_type_ == "WATER") {
+                                distr[pu.phase_pos[BlackoilPhases::Aqua]] = 1.0;
+                            } else if (wci_line.injector_type_ == "OIL") {
+                                distr[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
+                            } else if (wci_line.injector_type_ == "GAS") {
+                                distr[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                            } else {
+                                THROW("Injector type " << wci_line.injector_type_ << " not supported."
+                                      "WellsManager only supports WATER, OIL and GAS injector types.");
+                            }
                             ok = append_well_controls(RESERVOIR_RATE, wci_line.reservoir_flow_max_rate_,
                                                       distr, wix, w_);
                         }

tests/testFluid.DATA

@@ -0,0 +1,44 @@
+OIL
+WATER
+GAS
+
+FIELD
+
+PVTO
+--     Rs       Pbub       Bo        Vo
+      .0         14.7     1.0000     1.20  /
+      .165        400.    1.0120     1.17  /
+      .335        800.    1.0255     1.14  /
+      .500       1200.    1.0380     1.11  /
+      .665       1600.    1.0510     1.08  /
+      .828       2000.    1.0630     1.06  /
+      .985       2400.    1.0750     1.03  /
+     1.130       2800.    1.0870     1.00  /
+     1.270       3200.    1.0985      .98  /
+     1.390       3600.    1.1100      .95  /
+     1.500       4000.    1.1200      .94
+                 5000.    1.1189      .94  /
+ /
+
+PVDG
+--    Pg        Bg       Vg
+      14.7    178.08    .0125
+       400.     5.4777  .0130
+       800.     2.7392  .0135
+      1200.     1.8198  .0140
+      1600.     1.3648  .0145
+      2000.     1.0957  .0150
+      2400.     0.9099  .0155
+      2800.     0.7799  .0160
+      3200.     0.6871  .0165
+      3600.     0.6035  .0170
+      4000.     0.5432  .0175  /
+
+PVTW
+--Depth  Bw      Comp   Vw    Cv
+  3600. 1.0034  1.0E-6  0.96  0.0 /
+
+
+DENSITY
+-- Oil   Water   Gas
+ 44.98  63.01  0.0702  /