Some adjustments to equil initialisation.

- Saturations, phase pressures, and standard initialsation of RS and RV
   now agree to baseline.
 - Tables of RS and RV versus vertical depth (kw RSVD RVVD) have been
   hardcoded for testing (need new parser) and the calculations agree to
   baseline in the gas and oil zones.  In the water zone there is some
   differences: Our code computes saturated RS and RV using the final
   phase pressures (these are modified to be consistent with saturations
   and capillary pressures) while the baseline uses unmodified phase pressures.

opm/core/simulator/EquilibrationHelpers.hpp

@@ -175,7 +175,8 @@ namespace Opm
                  * depth and pressure @c press.
                  */
                 virtual double operator()(const double depth,
-                                          const double press) const = 0;
+                                          const double press,
+                                          const double sat = 0.0) const = 0;
             };
 
 
@@ -199,7 +200,8 @@ namespace Opm
                  */
                 double
                 operator()(const double /* depth */,
-                           const double /* press */) const
+                           const double /* press */,
+                           const double sat = 0.0) const
                 {
                     return 0.0;
                 }
@@ -216,14 +218,24 @@ namespace Opm
                 /**
                  * Constructor.
                  *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
                  * \param[in] depth Depth nodes.
                  * \param[in] rs Dissolved gas-oil ratio at @c depth.
                  */
-                RsVD(const std::vector<double>& depth,
+                RsVD(const BlackoilPropertiesInterface& props,
+                     const int cell,
+                     const std::vector<double>& depth,
                      const std::vector<double>& rs)
-                    : depth_(depth)
+                    : props_(props) 
+                    , cell_(cell)
+                    , depth_(depth)
                     , rs_(rs)
                 {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
                 }
 
                 /**
@@ -240,14 +252,111 @@ namespace Opm
                  */
                 double
                 operator()(const double depth,
-                           const double /* press */) const
+                           const double press,
+                           const double sat_gas = 0.0) const
                 {
-                    return linearInterpolation(depth_, rs_, depth);
+                    if (sat_gas > 0.0) {
+                        return satRs(press);
+                    } else {
+                        return std::min(satRs(press), linearInterpolation(depth_, rs_, depth));
+                    }
                 }
 
             private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
                 std::vector<double> depth_; /**< Depth nodes */
                 std::vector<double> rs_;    /**< Dissolved gas-oil ratio */
+                double z_[BlackoilPhases::MaxNumPhases];
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRs(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rs/Bo is in the gas row and oil column of A_.
+                    // 1/Bo is in the oil row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*opos + gpos] / A_[np*opos + opos];
+                }
+            };
+
+
+            /**
+             * Type that implements "vaporized oil-gas ratio"
+             * tabulated as a function of depth policy.  Data
+             * typically taken from keyword 'RVVD'.
+             */
+            class RvVD : public RsFunction {
+            public:
+                /**
+                 * Constructor.
+                 *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
+                 * \param[in] depth Depth nodes.
+                 * \param[in] rv Dissolved gas-oil ratio at @c depth.
+                 */
+                RvVD(const BlackoilPropertiesInterface& props,
+                     const int cell,
+                     const std::vector<double>& depth,
+                     const std::vector<double>& rv)
+                    : props_(props) 
+                    , cell_(cell)
+                    , depth_(depth)
+                    , rv_(rv)
+                {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
+                }
+
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RV
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RV
+                 * value.
+                 *
+                 * \return Vaporized oil-gas ratio (RV) at depth @c
+                 * depth and pressure @c press.
+                 */
+                double
+                operator()(const double depth,
+                           const double press,
+                           const double sat_oil = 0.0 ) const
+                {
+                    if (sat_oil > 0.0) {
+                        return satRv(press);
+                    } else {
+                        return std::min(satRv(press), linearInterpolation(depth_, rv_, depth));
+                    }
+                }
+
+            private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
+                std::vector<double> depth_; /**< Depth nodes */
+                std::vector<double> rv_;    /**< Vaporized oil-gas ratio */
+                double z_[BlackoilPhases::MaxNumPhases];
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRv(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rv/Bg is in the oil row and gas column of A_.
+                    // 1/Bg is in the gas row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*gpos + opos] / A_[np*gpos + gpos];
+                }
             };
 
 
@@ -298,9 +407,14 @@ namespace Opm
                  */
                 double
                 operator()(const double /* depth */,
-                           const double press) const
+                           const double press,
+                           const double sat_gas = 0.0) const
                 {                     
-                    return std::min(satRs(press), rs_sat_contact_);
+                    if (sat_gas > 0.0) {
+                        return satRs(press);
+                    } else {
+                        return std::min(satRs(press), rs_sat_contact_);
+                    }
                 }
 
             private:
@@ -323,6 +437,84 @@ namespace Opm
                 }
             };
 
+
+            /**
+             * Class that implements "vaporized oil-gas ratio" (Rv)
+             * as function of depth and pressure as follows:
+             *
+             *   1. The Rv at the gas-oil contact is equal to the
+             *      saturated Rv value, Rv_sat_contact.
+             *
+             *   2. The Rv elsewhere is equal to Rv_sat_contact, but
+             *      constrained to the saturated value as given by the
+             *      local pressure.
+             *
+             * This should yield Rv-values that are constant below the
+             * contact, and decreasing above the contact.
+             */
+            class RvSatAtContact : public RsFunction {
+            public:
+                /**
+                 * Constructor.
+                 *
+                 * \param[in] props      property object
+                 * \param[in] cell       any cell in the pvt region
+                 * \param[in] p_contact  oil pressure at the contact
+                 */
+                RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
+                    : props_(props), cell_(cell)
+                {
+                    auto pu = props_.phaseUsage();
+                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
+                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
+                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
+                    rv_sat_contact_ = satRv(p_contact);
+                }
+
+                /**
+                 * Function call.
+                 *
+                 * \param[in] depth Depth at which to calculate RV
+                 * value.
+                 *
+                 * \param[in] press Pressure at which to calculate RV
+                 * value.
+                 *
+                 * \return Dissolved oil-gas ratio (RV) at depth @c
+                 * depth and pressure @c press.
+                 */
+                double
+                operator()(const double /*depth*/,
+                           const double press,
+                           const double sat_oil = 0.0) const
+                {
+                    if (sat_oil > 0.0) {
+                        return satRv(press);
+                    } else {
+                        return std::min(satRv(press), rv_sat_contact_);
+                    }
+                }
+
+            private:
+                const BlackoilPropertiesInterface& props_;
+                const int cell_;
+                double z_[BlackoilPhases::MaxNumPhases];
+                double rv_sat_contact_;
+                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
+
+                double satRv(const double press) const
+                {
+                    props_.matrix(1, &press, z_, &cell_, A_, 0);
+                    // Rv/Bg is in the oil row and gas column of A_.
+                    // 1/Bg is in the gas row and column.
+                    // Recall also that it is stored in column-major order.
+                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
+                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
+                    const int np = props_.numPhases();
+                    return A_[np*gpos + opos] / A_[np*gpos + gpos];
+                }
+            };
+
         } // namespace Miscibility
 
 
@@ -355,6 +547,9 @@ namespace Opm
                 double depth;
                 double press;
             } main, woc, goc;
+            int live_oil_table_index;
+            int wet_gas_table_index;
+            int N;
         };
 
         /**

opm/core/simulator/initStateEquil.hpp

@@ -175,8 +175,8 @@ namespace Opm
         std::vector<double> computeRs(const UnstructuredGrid& grid,
                                       const CellRangeType& cells,
                                       const std::vector<double> oil_pressure,
-                                      const Miscibility::RsFunction& rs_func);
+                                      const Miscibility::RsFunction& rs_func,
-
+                                      const std::vector<double> gas_saturation);
 
         namespace DeckDependent {
 
@@ -205,8 +205,17 @@ namespace Opm
                                 ,
                                 { rec.gas_oil_contact_depth_   ,
                                   rec.gas_oil_cap_pressure_    }
+                                ,
+                                rec.live_oil_table_index_
+                                ,
+                                rec.wet_gas_table_index_
+                                ,
+                                rec.N_
                             };
-
+                        if (record.N != 0) {
+                            OPM_THROW(std::domain_error,
+                              "kw EQUIL, item 9: Only N=0 supported.");
+                        }
                         ret.push_back(record);
                     }
 
@@ -268,13 +277,24 @@ namespace Opm
                     // Create Rs functions.
                     rs_func_.reserve(rec.size());
                     if (deck.hasField("DISGAS")) {                    
-                        if (deck.hasField("RSVD")) {
+                        for (size_t i = 0; i < rec.size(); ++i) {
-                            // Rs has been specified as a function of depth.
+                            const int cell = *(eqlmap.cells(i + 1).begin());
-                            OPM_THROW(std::runtime_error, "Cannot initialise: RSVD field not read by EclipseGridParser class.");
+                            // TODO Check this ... 
-                        } else {
+                            //   The index <cell> is here picked as a representative for its equlibrium region,
-                            // Default initialisation: constant Rs below contact, saturated above.
+                            //   but is below used to identify PVT region.
-                            for (size_t i = 0; i < rec.size(); ++i) {
+                            //   This assumes that an eq-region has uniform pvt properties, is this always ok?
-                                const int cell = *(eqlmap.cells(i + 1).begin());
+                            //   For Norne this is trivial, as there is only one active pvt-region (PVTNUM=1 for all cells):                   
+                            if (rec[i].live_oil_table_index > 0) {
+                                if (deck.hasField("RSVD")) { 
+                                    // TODO When this kw is actually parsed, also check for proper number of available tables
+                                    // For now, just use dummy ...
+                                    std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
+                                    std::vector<double> rs; rs.push_back(0.0); rs.push_back(100.0);
+                                    rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, depth, rs));
+                                } else {
+                                    OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
+                                }
+                            } else {
                                 if (rec[i].goc.depth != rec[i].main.depth) {
                                     OPM_THROW(std::runtime_error,
                                               "Cannot initialise: when no explicit RSVD table is given, \n"
@@ -291,6 +311,38 @@ namespace Opm
                         }
                     }                    
 
+                    rv_func_.reserve(rec.size());
+                    if (deck.hasField("VAPOIL")) {                    
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            const int cell = *(eqlmap.cells(i + 1).begin());
+                            // TODO Similar as above: eq-region vs pvt-region ...                   
+                            if (rec[i].wet_gas_table_index > 0) {
+                                if (deck.hasField("RVVD")) { 
+                                    // TODO When this kw is actually parsed, also check for proper number of available tables
+                                    // For now, just use dummy ...
+                                    std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
+                                    std::vector<double> rv; rv.push_back(0.0); rv.push_back(0.0001);
+                                    rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, depth, rv));
+                                } else {
+                                    OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
+                                }
+                            } else {
+                                if (rec[i].goc.depth != rec[i].main.depth) {
+                                    OPM_THROW(std::runtime_error,
+                                              "Cannot initialise: when no explicit RVVD table is given, \n"
+                                              "datum depth must be at the gas-oil-contact. "
+                                              "In EQUIL region " << (i + 1) << "  (counting from 1), this does not hold.");
+                                }
+                                const double p_contact = rec[i].main.press + rec[i].goc.press;
+                                rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
+                            }
+                        }
+                    } else {
+                        for (size_t i = 0; i < rec.size(); ++i) {
+                            rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
+                        }
+                    }
+                    
                     // Compute pressures, saturations, rs and rv factors.
                     calcPressSatRsRv(eqlmap, rec, props, G, grav);
 
@@ -311,6 +363,7 @@ namespace Opm
                 typedef EquilReg<RhoCalc> EqReg;
 
                 std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
+                std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
 
                 PVec pp_;
                 PVec sat_;
@@ -335,13 +388,14 @@ namespace Opm
 
                         const int repcell = *cells.begin();
                         const RhoCalc calc(props, repcell);
-
                         const EqReg eqreg(rec[r], calc,
-                                          rs_func_[r], std::make_shared<NoMix>(),
+                                          rs_func_[r], rv_func_[r],
                                           props.phaseUsage());
                                           
-                        const PVec press = phasePressures(G, eqreg, cells, grav);
+                        PVec press = phasePressures(G, eqreg, cells, grav);
+                        
                         const PVec sat = phaseSaturations(eqreg, cells, props, press);
+                        
                         const int np = props.numPhases();
                         for (int p = 0; p < np; ++p) {
                             copyFromRegion(press[p], cells, pp_[p]);
@@ -350,8 +404,9 @@ namespace Opm
                         if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
                             && props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
                             const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
-                            const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]));
+                            const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
-                            const Vec rv(cells.size(), 0.0);
+                            const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
+                            const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
                             copyFromRegion(rs, cells, rs_);
                             copyFromRegion(rv, cells, rv_);
                         }

opm/core/simulator/initStateEquil_impl.hpp

@@ -580,7 +580,7 @@ namespace Opm
         phaseSaturations(const Region&           reg,
                          const CellRange&        cells,
                          const BlackoilPropertiesInterface& props,
-                         const std::vector< std::vector<double> >& phase_pressures)
+                         std::vector< std::vector<double> >& phase_pressures)
         {
             const double z0   = reg.datum();
             const double zwoc = reg.zwoc ();
@@ -621,6 +621,7 @@ namespace Opm
                     sg = satFromPc(props, gaspos, cell, pcog, increasing);
                     phase_saturations[gaspos][local_index] = sg;
                 }
+                bool overlap = false;
                 if (gas && water && (sg + sw > 1.0)) {
                     // Overlapping gas-oil and oil-water transition
                     // zones can lead to unphysical saturations when
@@ -631,8 +632,32 @@ namespace Opm
                     sg = 1.0 - sw;
                     phase_saturations[waterpos][local_index] = sw;
                     phase_saturations[gaspos][local_index] = sg;
+                    overlap = true;
                 }
                 phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
+                
+                // Adjust phase pressures for max and min saturation ...
+                double pc[BlackoilPhases::MaxNumPhases];
+                double sat[BlackoilPhases::MaxNumPhases];
+                if (sw > smax[waterpos]-1.0e-6) {
+                    sat[waterpos] = smax[waterpos];
+                    props.capPress(1, sat, &cell, pc, 0);                   
+                    phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
+                } else if (overlap || sg > smax[gaspos]-1.0e-6) {
+                    sat[gaspos] = smax[gaspos];
+                    props.capPress(1, sat, &cell, pc, 0);                   
+                    phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
+                }
+                if (sg < smin[gaspos]+1.0e-6) {
+                    sat[gaspos] = smin[gaspos];
+                    props.capPress(1, sat, &cell, pc, 0);
+                    phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
+                }
+                if (sw < smin[waterpos]+1.0e-6) {
+                    sat[waterpos] = smin[waterpos];
+                    props.capPress(1, sat, &cell, pc, 0);
+                    phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
+                }
             }
             return phase_saturations;
         }
@@ -659,20 +684,19 @@ namespace Opm
         std::vector<double> computeRs(const UnstructuredGrid& grid,
                                       const CellRangeType& cells,
                                       const std::vector<double> oil_pressure,
-                                      const Miscibility::RsFunction& rs_func)
+                                      const Miscibility::RsFunction& rs_func,
+                                      const std::vector<double> gas_saturation)
         {
             assert(grid.dimensions == 3);
             std::vector<double> rs(cells.size());
             int count = 0;
             for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
                 const double depth = grid.cell_centroids[3*(*it) + 2];
-                rs[count] = rs_func(depth, oil_pressure[count]);
+                rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]);
             }
             return rs;
         }
 
-
-
     } // namespace Equil