diff --git a/examples/compute_initial_state.cpp b/examples/compute_initial_state.cpp index 387c1fef..3b118c15 100644 --- a/examples/compute_initial_state.cpp +++ b/examples/compute_initial_state.cpp @@ -95,7 +95,7 @@ try // Initialisation. BlackoilState state; - initStateEquil(grid, props, deck, grav, state); + initStateEquil(grid, props, deck, eclipseState, grav, state); // Output. const std::string output_dir = param.getDefault("output_dir", "output"); diff --git a/opm/core/props/BlackoilPropertiesBasic.cpp b/opm/core/props/BlackoilPropertiesBasic.cpp index ab9b238b..9634d9aa 100644 --- a/opm/core/props/BlackoilPropertiesBasic.cpp +++ b/opm/core/props/BlackoilPropertiesBasic.cpp @@ -239,5 +239,16 @@ namespace Opm } + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + void BlackoilPropertiesBasic::swatInitScaling(const int cell, + const double pcow, + double & swat) + { + OPM_THROW(std::runtime_error, "BlackoilPropertiesBasic::swatInitScaling() -- not implemented."); + } + } // namespace Opm diff --git a/opm/core/props/BlackoilPropertiesBasic.hpp b/opm/core/props/BlackoilPropertiesBasic.hpp index b85991ed..5c0472b9 100644 --- a/opm/core/props/BlackoilPropertiesBasic.hpp +++ b/opm/core/props/BlackoilPropertiesBasic.hpp @@ -163,9 +163,9 @@ namespace Opm double* dpcds) const; -/// Obtain the range of allowable saturation values. - /// In cell cells[i], saturation of phase p is allowed to be - /// in the interval [smin[i*P + p], smax[i*P + p]]. + /// Obtain the range of allowable saturation values. + /// In cell cells[i], saturation of phase p is allowed to be + /// in the interval [smin[i*P + p], smax[i*P + p]]. /// \param[in] n Number of data points. /// \param[in] cells Array of n cell indices. /// \param[out] smin Array of nP minimum s values, array must be valid before calling. @@ -175,6 +175,16 @@ namespace Opm double* smin, double* smax) const; + + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + virtual void swatInitScaling(const int cell, + const double pcow, + double & swat); + + private: RockBasic rock_; PvtPropertiesBasic pvt_; diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index 9251ecdc..16cec166 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -303,7 +303,18 @@ namespace Opm { satprops_->satRange(n, cells, smin, smax); } - + + + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + void BlackoilPropertiesFromDeck::swatInitScaling(const int cell, + const double pcow, + double & swat) + { + satprops_->swatInitScaling(cell, pcow, swat); + } } // namespace Opm diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 069405c4..3c71bda3 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -215,6 +215,15 @@ namespace Opm const int* cells, double* smin, double* smax) const; + + + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + virtual void swatInitScaling(const int cell, + const double pcow, + double & swat); private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const diff --git a/opm/core/props/BlackoilPropertiesInterface.hpp b/opm/core/props/BlackoilPropertiesInterface.hpp index 08edc8ae..8a25e344 100644 --- a/opm/core/props/BlackoilPropertiesInterface.hpp +++ b/opm/core/props/BlackoilPropertiesInterface.hpp @@ -160,6 +160,16 @@ namespace Opm const int* cells, double* smin, double* smax) const = 0; + + + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + virtual void swatInitScaling(const int cell, + const double pcow, + double & swat) = 0; + }; diff --git a/opm/core/props/satfunc/SatFuncBase.hpp b/opm/core/props/satfunc/SatFuncBase.hpp index dd80a10b..cb1cab3c 100644 --- a/opm/core/props/satfunc/SatFuncBase.hpp +++ b/opm/core/props/satfunc/SatFuncBase.hpp @@ -55,6 +55,7 @@ namespace Opm double krSlopeCrit; double scaleKr(double s, double kr, double krsr_tab) const; double scaleKrDeriv(double s, double krDeriv) const; // Returns scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s)) + double pcFactor; // Scaling factor for capillary pressure. void printMe(std::ostream & out); }; @@ -99,6 +100,8 @@ namespace Opm double sogcr_; // Critical oil-in-gas-and-connate-water saturation. double krorw_; // Oil relperm at critical water saturation. double krorg_; // Oil relperm at critical gas saturation. + double pcwmax_; // Max oil-water capillary pressure. + double pcgmax_; // Max gas-oil capillary pressure. protected: PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. @@ -182,6 +185,7 @@ namespace Opm break; } } + pcwmax_ = pcow.front(); } if (phase_usage.phase_used[Vapour]) { Opm::SgofTable sgof(deck->getKeyword("SGOF"), table_num); @@ -234,6 +238,7 @@ namespace Opm break; } } + pcgmax_ = pcog.back(); } diff --git a/opm/core/props/satfunc/SatFuncGwseg.hpp b/opm/core/props/satfunc/SatFuncGwseg.hpp index a8dfb33f..aaec8fa9 100644 --- a/opm/core/props/satfunc/SatFuncGwseg.hpp +++ b/opm/core/props/satfunc/SatFuncGwseg.hpp @@ -553,12 +553,12 @@ namespace Opm if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; double _sw = epst->wat.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); - pc[pos] = this->pcow_(_sw); + pc[pos] = epst->wat.pcFactor*this->pcow_(_sw); } if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); - pc[pos] = this->pcog_(_sg); + pc[pos] = epst->gas.pcFactor*this->pcog_(_sg); } } @@ -602,16 +602,16 @@ namespace Opm if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; double _sw = epst->wat.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); - pc[pos] = this->pcow_(s[pos]); + pc[pos] = epst->wat.pcFactor*this->pcow_(s[pos]); double _dsdsw = epst->wat.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]); - dpcds[np*pos + pos] = _dsdsw*this->pcow_.derivative(_sw); + dpcds[np*pos + pos] = epst->wat.pcFactor*_dsdsw*this->pcow_.derivative(_sw); } if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]); - pc[pos] = this->pcog_(_sg); + pc[pos] = epst->gas.pcFactor*this->pcog_(_sg); double _dsdsg = epst->gas.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]); - dpcds[np*pos + pos] = _dsdsg*this->pcog_.derivative(_sg); + dpcds[np*pos + pos] = epst->gas.pcFactor*_dsdsg*this->pcog_.derivative(_sg); } } diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index d9c96ea3..8db30ac3 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -124,7 +124,7 @@ namespace Opm const int* cells, double* smin, double* smax) const; - + /// Update saturation state for the hysteresis tracking /// \param[in] n Number of data points. /// \param[in] s Array of nP saturation values. @@ -132,6 +132,14 @@ namespace Opm const int* cells, const double* s); + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + void swatInitScaling(const int cell, + const double pcow, + double & swat); + private: PhaseUsage phase_usage_; std::vector satfuncset_; @@ -178,13 +186,15 @@ namespace Opm const double s0_tab, const double krsr_tab, const double krmax_tab, + const double pcmax_tab, const std::vector& sl, const std::vector& scr, const std::vector& su, const std::vector& sxcr, const std::vector& s0, const std::vector& krsr, - const std::vector& krmax); + const std::vector& krmax, + const std::vector& pcmax); bool columnIsMasked_(Opm::DeckConstPtr deck, const std::string& keywordName, diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index b1124ffa..f51dde88 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -412,6 +412,39 @@ namespace Opm } } + + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + template + void SaturationPropsFromDeck::swatInitScaling(const int cell, + const double pcow, + double & swat) + { + if (phase_usage_.phase_used[Aqua]) { + // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74 + if (swat <= eps_transf_[cell].wat.smin) { + swat = eps_transf_[cell].wat.smin; + } else if (pcow < 1.0e-8) { + swat = eps_transf_[cell].wat.smax; + } else { + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int np = phase_usage_.num_phases; + double s[np]; + s[wpos] = swat; + double pc[np]; + funcForCell(cell).evalPc(s, pc, &(eps_transf_[cell])); + if (pc[wpos] > 1.0e-8) { + eps_transf_[cell].wat.pcFactor *= pcow/pc[wpos]; + } + } + } else { + OPM_THROW(std::runtime_error, "swatInitScaling: no water phase! "); + } + } + + // Map the cell number to the correct function set. template const typename SaturationPropsFromDeck::Funcs& @@ -431,6 +464,8 @@ namespace Opm { std::vector swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr; std::vector krw, krg, kro, krwr, krgr, krorw, krorg; + std::vector pcw, pcg; + const std::vector dummy; // Initialize saturation scaling parameter initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, std::string("SWL"), swl); @@ -462,6 +497,10 @@ namespace Opm std::string("KRORW"), krorw); initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, std::string("KRORG"), krorg); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("PCW"), pcw); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("PCG"), pcg); eps_transf_.resize(number_of_cells); @@ -474,31 +513,95 @@ namespace Opm for (int cell = 0; cell < number_of_cells; ++cell) { if (oilWater) { // ### krw - initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], - funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw); + initEPSParam(cell, eps_transf_[cell].wat, false, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, + -1.0, + funcForCell(cell).krwr_, + funcForCell(cell).krwmax_, + funcForCell(cell).pcwmax_, + swl, swcr, swu, sowcr, sgl, krwr, krw, pcw); // ### krow - initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], - funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro); + initEPSParam(cell, eps_transf_[cell].watoil, true, + 0.0, + funcForCell(cell).sowcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + -1.0, + funcForCell(cell).krorw_, + funcForCell(cell).kromax_, + 0.0, + swl, sowcr, swl, swcr, sgl, krorw, kro, dummy); } else if (oilGas) { // ### krg - initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], - funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg); + initEPSParam(cell, eps_transf_[cell].gas, false, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, + -1.0, + funcForCell(cell).krgr_, + funcForCell(cell).krgmax_, + funcForCell(cell).pcgmax_, + sgl, sgcr, sgu, sogcr, swl, krgr, krg, pcg); // ### krog - initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], - funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro); + initEPSParam(cell, eps_transf_[cell].gasoil, true, + 0.0, + funcForCell(cell).sogcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + -1.0, + funcForCell(cell).krorg_, + funcForCell(cell).kromax_, + 0.0, + sgl, sogcr, sgl, sgcr, swl, krorg, kro, dummy); } else if (threephase) { // ### krw - initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_, - funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw); - // ### krow - initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, - funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro); - // ### krg - initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_, - funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg); - // ### krog - initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, - funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro); + initEPSParam(cell, eps_transf_[cell].wat, false, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).krwr_, + funcForCell(cell).krwmax_, + funcForCell(cell).pcwmax_, + swl, swcr, swu, sowcr, sgl, krwr, krw, pcw); + // ### krow + initEPSParam(cell, eps_transf_[cell].watoil, true, + 0.0, + funcForCell(cell).sowcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).krorw_, + funcForCell(cell).kromax_, + 0.0, + swl, sowcr, swl, swcr, sgl, krorw, kro, dummy); + // ### krg + initEPSParam(cell, eps_transf_[cell].gas, false, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).krgr_, + funcForCell(cell).krgmax_, + funcForCell(cell).pcgmax_, + sgl, sgcr, sgu, sogcr, swl, krgr, krg, pcg); + // ### krog + initEPSParam(cell, eps_transf_[cell].gasoil, true, + 0.0, + funcForCell(cell).sogcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).krorg_, + funcForCell(cell).kromax_, + 0.0, + sgl, sogcr, sgl, sgcr, swl, krorg, kro, dummy); } } } @@ -514,6 +617,8 @@ namespace Opm { std::vector iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr; std::vector ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg; + std::vector ipcw, ipcg; + const std::vector dummy; // Initialize hysteresis saturation scaling parameters initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, std::string("ISWL"), iswl); @@ -545,6 +650,10 @@ namespace Opm std::string("IKRORW"), ikrorw); initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, std::string("IKRORG"), ikrorg); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IPCW"), ipcw); + initEPSKey(deck, number_of_cells, global_cell, begin_cell_centroid, dimensions, + std::string("IPCG"), ipcg); eps_transf_hyst_.resize(number_of_cells); sat_hyst_.resize(number_of_cells); @@ -557,32 +666,96 @@ namespace Opm for (int cell = 0; cell < number_of_cells; ++cell) { if (oilWater) { - // ### krw - initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], - funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw); - // ### krow - initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], - funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro); - } else if (oilGas) { - // ### krg - initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], - funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg); - // ### krog - initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], - funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro); - } else if (threephase) { - // ### krw - initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_, - funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw); - // ### krow - initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, - funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro); - // ### krg - initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_, - funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg); - // ### krog - initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, - funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro); + // ### krw + initEPSParam(cell, eps_transf_hyst_[cell].wat, false, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, + -1.0, + funcForCell(cell).krwr_, + funcForCell(cell).krwmax_, + funcForCell(cell).pcwmax_, + iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw, ipcw); + // ### krow + initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, + 0.0, + funcForCell(cell).sowcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + -1.0, + funcForCell(cell).krorw_, + funcForCell(cell).kromax_, + 0.0, + iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro, dummy); + } else if (oilGas) { + // ### krg + initEPSParam(cell, eps_transf_hyst_[cell].gas, false, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, + -1.0, + funcForCell(cell).krgr_, + funcForCell(cell).krgmax_, + funcForCell(cell).pcgmax_, + isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg, ipcg); + // ### krog + initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, + 0.0, + funcForCell(cell).sogcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + -1.0, + funcForCell(cell).krorg_, + funcForCell(cell).kromax_, + 0.0, + isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro, dummy); + } else if (threephase) { + // ### krw + initEPSParam(cell, eps_transf_hyst_[cell].wat, false, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).krwr_, + funcForCell(cell).krwmax_, + funcForCell(cell).pcwmax_, + iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw, ipcw); + // ### krow + initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, + 0.0, + funcForCell(cell).sowcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).krorw_, + funcForCell(cell).kromax_, + 0.0, + iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro, dummy); + // ### krg + initEPSParam(cell, eps_transf_hyst_[cell].gas, false, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).krgr_, + funcForCell(cell).krgmax_, + funcForCell(cell).pcgmax_, + isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg, ipcg); + // ### krog + initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, + 0.0, + funcForCell(cell).sogcr_, + funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, + funcForCell(cell).smin_[wpos], + funcForCell(cell).krorg_, + funcForCell(cell).kromax_, + 0.0, + isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro, dummy); } } } @@ -747,6 +920,16 @@ namespace Opm param_col[table_num] = enkrvd.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg } } + } else if (useKeyword && (keyword[0] == 'P' || keyword[1] == 'P') ) { + if (useAqua && (keyword == std::string("PCW") || keyword == std::string("IPCW")) ) { + scaleparam.resize(number_of_cells); + for (int i=0; i no scaling) const std::vector& sl, // For krow/krog calculations this is not used const std::vector& scr, const std::vector& su, // For krow/krog calculations this is SWL/SGL const std::vector& sxcr, const std::vector& s0, const std::vector& krsr, - const std::vector& krmax) + const std::vector& krmax, + const std::vector& pcmax) // For krow/krog calculations this is not used { if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) { data.doNotScale = true; @@ -876,6 +1061,12 @@ namespace Opm } } + if (std::fabs(pcmax_tab) < 1.0e-8 || pcmax.empty() || pcmax_tab*pcmax[cell] < 0.0) { + data.pcFactor = 1.0; + } else { + data.pcFactor = pcmax[cell]/pcmax_tab; + } + } } // namespace Opm diff --git a/opm/core/props/satfunc/SaturationPropsInterface.hpp b/opm/core/props/satfunc/SaturationPropsInterface.hpp index 04f6730e..2b54d39e 100644 --- a/opm/core/props/satfunc/SaturationPropsInterface.hpp +++ b/opm/core/props/satfunc/SaturationPropsInterface.hpp @@ -80,7 +80,14 @@ namespace Opm virtual void updateSatHyst(const int n, const int* cells, const double* s) = 0; - + + /// Update capillary pressure scaling according to pressure diff. and initial water saturation. + /// \param[in] cell Cell index. + /// \param[in] pcow P_oil - P_water. + /// \param[in/out] swat Water saturation. / Possibly modified Water saturation. + virtual void swatInitScaling(const int cell, + const double pcow, + double & swat) = 0; }; diff --git a/opm/core/simulator/initStateEquil.hpp b/opm/core/simulator/initStateEquil.hpp index a4ec1ee4..b869aa97 100644 --- a/opm/core/simulator/initStateEquil.hpp +++ b/opm/core/simulator/initStateEquil.hpp @@ -28,6 +28,7 @@ #include #include #include +#include #include #include @@ -62,6 +63,7 @@ namespace Opm void initStateEquil(const UnstructuredGrid& grid, const BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, + const Opm::EclipseStateConstPtr eclipseState, const double gravity, BlackoilState& state); @@ -150,8 +152,9 @@ namespace Opm std::vector< std::vector > phaseSaturations(const Region& reg, const CellRange& cells, - const BlackoilPropertiesInterface& props, - const std::vector< std::vector >& phase_pressures); + BlackoilPropertiesInterface& props, + const std::vector swat_init, + std::vector< std::vector >& phase_pressures); @@ -229,13 +232,14 @@ namespace Opm inline std::vector equilnum(const Opm::DeckConstPtr deck, + const Opm::EclipseStateConstPtr eclipseState, const UnstructuredGrid& G ) { std::vector eqlnum; if (deck->hasKeyword("EQLNUM")) { eqlnum.resize(G.number_of_cells); const std::vector& e = - deck->getKeyword("EQLNUM")->getIntData(); + eclipseState->getIntGridProperty("EQLNUM")->getData(); const int* gc = G.global_cell; for (int cell = 0; cell < G.number_of_cells; ++cell) { const int deck_pos = (gc == NULL) ? cell : gc[cell]; @@ -254,8 +258,9 @@ namespace Opm class InitialStateComputer { public: - InitialStateComputer(const BlackoilPropertiesInterface& props, + InitialStateComputer(BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, + const Opm::EclipseStateConstPtr eclipseState, const UnstructuredGrid& G , const double grav = unit::gravity) : pp_(props.numPhases(), @@ -269,7 +274,7 @@ namespace Opm const std::vector rec = getEquil(deck); // Create (inverse) region mapping. - const RegionMapping<> eqlmap(equilnum(deck, G)); + const RegionMapping<> eqlmap(equilnum(deck, eclipseState, G)); // Create Rs functions. rs_func_.reserve(rec.size()); @@ -333,6 +338,18 @@ namespace Opm rv_func_.push_back(std::make_shared()); } } + + + // Check for presence of kw SWATINIT + if (deck->hasKeyword("SWATINIT")) { + const std::vector& swat_init = eclipseState->getDoubleGridProperty("SWATINIT")->getData(); + swat_init_.resize(G.number_of_cells); + const int* gc = G.global_cell; + for (int c = 0; c < G.number_of_cells; ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + swat_init_[c] = swat_init[deck_pos]; + } + } // Compute pressures, saturations, rs and rv factors. calcPressSatRsRv(eqlmap, rec, props, G, grav); @@ -360,15 +377,18 @@ namespace Opm PVec sat_; Vec rs_; Vec rv_; + Vec swat_init_; template void - calcPressSatRsRv(const RMap& reg , - const std::vector< EquilRecord >& rec , - const Opm::BlackoilPropertiesInterface& props, - const UnstructuredGrid& G , + calcPressSatRsRv(const RMap& reg , + const std::vector< EquilRecord >& rec , + Opm::BlackoilPropertiesInterface& props, + const UnstructuredGrid& G , const double grav) { + typedef Miscibility::NoMixing NoMix; + for (typename RMap::RegionId r = 0, nr = reg.numRegions(); r < nr; ++r) @@ -383,7 +403,7 @@ namespace Opm PVec press = phasePressures(G, eqreg, cells, grav); - const PVec sat = phaseSaturations(eqreg, cells, props, press); + const PVec sat = phaseSaturations(eqreg, cells, props, swat_init_, press); const int np = props.numPhases(); for (int p = 0; p < np; ++p) { diff --git a/opm/core/simulator/initStateEquil_impl.hpp b/opm/core/simulator/initStateEquil_impl.hpp index ba2df1ed..611c5cd4 100644 --- a/opm/core/simulator/initStateEquil_impl.hpp +++ b/opm/core/simulator/initStateEquil_impl.hpp @@ -22,6 +22,7 @@ #include #include +#include #include #include @@ -579,7 +580,8 @@ namespace Opm std::vector< std::vector > phaseSaturations(const Region& reg, const CellRange& cells, - const BlackoilPropertiesInterface& props, + BlackoilPropertiesInterface& props, + const std::vector swat_init, std::vector< std::vector >& phase_pressures) { const double z0 = reg.datum(); @@ -610,8 +612,14 @@ namespace Opm double sw = 0.0; if (water) { const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; - sw = satFromPc(props, waterpos, cell, pcov); - phase_saturations[waterpos][local_index] = sw; + if (swat_init.empty()) { // Invert Pc to find sw + sw = satFromPc(props, waterpos, cell, pcov); + phase_saturations[waterpos][local_index] = sw; + } else { // Scale Pc to reflect imposed sw + sw = swat_init[cell]; + props.swatInitScaling(cell, pcov, sw); + phase_saturations[waterpos][local_index] = sw; + } } double sg = 0.0; if (gas) { @@ -621,18 +629,28 @@ 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 // treated as above. Must recalculate using gas-water // capillary pressure. const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index]; + if (! swat_init.empty()) { + // Re-scale Pc to reflect imposed sw for vanishing oil phase. + // This seems consistent with ecl, and fails to honour + // swat_init in case of non-trivial gas-oil cap pressure. + props.swatInitScaling(cell, pcgw, sw); + } sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw); sg = 1.0 - sw; phase_saturations[waterpos][local_index] = sw; phase_saturations[gaspos][local_index] = sg; - overlap = true; + // Adjust oil pressure according to gas saturation and cap pressure + double pc[BlackoilPhases::MaxNumPhases]; + double sat[BlackoilPhases::MaxNumPhases]; + sat[gaspos] = sg; + props.capPress(1, sat, &cell, pc, 0); + phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos]; } phase_saturations[oilpos][local_index] = 1.0 - sw - sg; @@ -643,7 +661,7 @@ namespace Opm 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) { + } else if (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]; @@ -736,13 +754,14 @@ namespace Opm * \param[in] gravity Acceleration of gravity, assumed to be in Z direction. */ void initStateEquil(const UnstructuredGrid& grid, - const BlackoilPropertiesInterface& props, + BlackoilPropertiesInterface& props, const Opm::DeckConstPtr deck, + const Opm::EclipseStateConstPtr eclipseState, const double gravity, BlackoilState& state) { typedef Equil::DeckDependent::InitialStateComputer ISC; - ISC isc(props, deck, grid, gravity); + ISC isc(props, deck, eclipseState, grid, gravity); const auto pu = props.phaseUsage(); const int ref_phase = pu.phase_used[BlackoilPhases::Liquid] ? pu.phase_pos[BlackoilPhases::Liquid] @@ -751,7 +770,7 @@ namespace Opm state.saturation() = convertSats(isc.saturation()); state.gasoilratio() = isc.rs(); state.rv() = isc.rv(); - // TODO: state.surfacevol() must be computed from s, rs, rv. + initBlackoilSurfvolUsingRSorRV(grid, props, state); } diff --git a/tests/test_equil.cpp b/tests/test_equil.cpp index 7084b80f..83947855 100644 --- a/tests/test_equil.cpp +++ b/tests/test_equil.cpp @@ -338,7 +338,7 @@ BOOST_AUTO_TEST_CASE (DeckAllDead) Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA"); Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, *grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, *grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, *grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells); @@ -419,7 +419,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillary) Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 10.0); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 10.0); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -459,7 +459,7 @@ BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap) Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -521,7 +521,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveOil) Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -600,7 +600,7 @@ BOOST_AUTO_TEST_CASE (DeckWithLiveGas) Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells); @@ -682,7 +682,7 @@ BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD) Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, false); - Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, grid, 9.80665); + Opm::Equil::DeckDependent::InitialStateComputer comp(props, deck, eclipseState, grid, 9.80665); const auto& pressures = comp.press(); BOOST_REQUIRE(pressures.size() == 3); BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);