From 14e271f0496cba7ff39dca5d2edeadbfa0ff0660 Mon Sep 17 00:00:00 2001 From: osae Date: Wed, 26 Mar 2014 14:08:39 +0100 Subject: [PATCH] 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 | 213 +++++++++++++++++++- opm/core/simulator/initStateEquil.hpp | 91 +++++++-- opm/core/simulator/initStateEquil_impl.hpp | 34 +++- 3 files changed, 306 insertions(+), 32 deletions(-) diff --git a/opm/core/simulator/EquilibrationHelpers.hpp b/opm/core/simulator/EquilibrationHelpers.hpp index 0d4d0819..d89dd8c3 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 3767563f..566f2b00 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 c075f4ec..6c991215 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