diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 0d4d08196..d89dd8c3c 100644 --- a/opm/core/simulator/EquilibrationHelpers.hpp +++ b/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& depth, + RsVD(const BlackoilPropertiesInterface& props, + const int cell, + const std::vector& depth, const std::vector& 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 depth_; /**< Depth nodes */ std::vector 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& depth, + const std::vector& 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 depth_; /**< Depth nodes */ + std::vector 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 - { - return std::min(satRs(press), rs_sat_contact_); + const double press, + const double sat_gas = 0.0) const + { + 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; }; /** diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index 3767563f3..566f2b00f 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -175,8 +175,8 @@ namespace Opm std::vector computeRs(const UnstructuredGrid& grid, const CellRangeType& cells, const std::vector oil_pressure, - const Miscibility::RsFunction& rs_func); - + const Miscibility::RsFunction& rs_func, + const std::vector 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); } @@ -267,14 +276,25 @@ namespace Opm // Create Rs functions. rs_func_.reserve(rec.size()); - if (deck.hasField("DISGAS")) { - if (deck.hasField("RSVD")) { - // Rs has been specified as a function of depth. - OPM_THROW(std::runtime_error, "Cannot initialise: RSVD field not read by EclipseGridParser class."); - } else { - // Default initialisation: constant Rs below contact, saturated above. - for (size_t i = 0; i < rec.size(); ++i) { - const int cell = *(eqlmap.cells(i + 1).begin()); + if (deck.hasField("DISGAS")) { + for (size_t i = 0; i < rec.size(); ++i) { + const int cell = *(eqlmap.cells(i + 1).begin()); + // TODO Check this ... + // The index is here picked as a representative for its equlibrium region, + // but is below used to identify PVT region. + // This assumes that an eq-region has uniform pvt properties, is this always ok? + // 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 depth; depth.push_back(0.0); depth.push_back(100.0); + std::vector rs; rs.push_back(0.0); rs.push_back(100.0); + rs_func_.push_back(std::make_shared(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" @@ -289,8 +309,40 @@ namespace Opm for (size_t i = 0; i < rec.size(); ++i) { rs_func_.push_back(std::make_shared()); } - } + } + 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 depth; depth.push_back(0.0); depth.push_back(100.0); + std::vector rv; rv.push_back(0.0); rv.push_back(0.0001); + rv_func_.push_back(std::make_shared(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(props, cell, p_contact)); + } + } + } else { + for (size_t i = 0; i < rec.size(); ++i) { + rv_func_.push_back(std::make_shared()); + } + } + // Compute pressures, saturations, rs and rv factors. calcPressSatRsRv(eqlmap, rec, props, G, grav); @@ -311,6 +363,7 @@ namespace Opm typedef EquilReg EqReg; std::vector< std::shared_ptr > rs_func_; + std::vector< std::shared_ptr > 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(), + 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 Vec rv(cells.size(), 0.0); + const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour]; + 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_); } diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index c075f4ecb..6c9912151 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/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 >& phase_pressures) + std::vector< std::vector >& 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 computeRs(const UnstructuredGrid& grid, const CellRangeType& cells, const std::vector oil_pressure, - const Miscibility::RsFunction& rs_func) + const Miscibility::RsFunction& rs_func, + const std::vector gas_saturation) { assert(grid.dimensions == 3); std::vector 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