From fa470063ca28f9d2e3b323aafa492219dd3fc7f5 Mon Sep 17 00:00:00 2001 From: osae Date: Tue, 28 Jan 2014 16:36:55 +0100 Subject: [PATCH] Saturation table scaling and hysteresis. This commit implements some additional scaling keywords. This includes the ISWL-family that provide hysteresis behaviour via alternative scaling of the tables. The old parser has been somewhat extended for testing purposes. The commit also includes a slight refactoring of the SatFunc-family where a new base class has been introduced. --- opm/core/io/eclipse/EclipseGridParser.cpp | 32 +- opm/core/io/eclipse/SpecialEclipseFields.hpp | 39 +- opm/core/props/satfunc/SatFuncBase.hpp | 312 +++++++++ opm/core/props/satfunc/SatFuncGwseg.cpp | 503 +------------- opm/core/props/satfunc/SatFuncGwseg.hpp | 643 +++++++++++++++-- opm/core/props/satfunc/SatFuncSimple.cpp | 646 +++++------------- opm/core/props/satfunc/SatFuncSimple.hpp | 306 +++++++-- opm/core/props/satfunc/SatFuncStone2.cpp | 389 +---------- opm/core/props/satfunc/SatFuncStone2.hpp | 216 ++++-- .../props/satfunc/SaturationPropsFromDeck.hpp | 50 +- .../satfunc/SaturationPropsFromDeck_impl.hpp | 581 ++++++++++------ 11 files changed, 1900 insertions(+), 1817 deletions(-) create mode 100644 opm/core/props/satfunc/SatFuncBase.hpp diff --git a/opm/core/io/eclipse/EclipseGridParser.cpp b/opm/core/io/eclipse/EclipseGridParser.cpp index f0106cc3..d40db19c 100644 --- a/opm/core/io/eclipse/EclipseGridParser.cpp +++ b/opm/core/io/eclipse/EclipseGridParser.cpp @@ -97,13 +97,19 @@ namespace EclipseKeywords string("SHEARMOD"), string("POISSONMOD"), string("PWAVEMOD"), string("MULTPV"), string("PRESSURE"), string("SGAS"), string("SWAT"), string("SOIL"), string("RS"), - string("RV"), string("DXV"), string("DYV"), string("DZV"), string("DEPTHZ"), string("TOPS"), string("MAPAXES"), - string("SWCR"), string("SWL"), string("SWU"), - string("SOWCR"), string("KRW"), string("KRWR"), - string("KRO"), string("KRORW"), string("NTG"), - string("RHO") + string("SWL"), string("SWCR"), string("SWU"), + string("SGL"), string("SGCR"), string("SGU"), + string("SOWCR"), string("SOGCR"), string("KRW"), + string("KRWR"), string("KRG"), string("KRGR"), + string("KRO"), string("KRORW"), string("KRORG"), + string("ISWL"), string("ISWCR"), string("ISWU"), + string("ISGL"), string("ISGCR"), string("ISGU"), + string("ISOWCR"), string("ISOGCR"), string("IKRW"), + string("IKRWR"), string("IKRG"), string("IKRGR"), + string("IKRO"), string("IKRORW"), string("IKRORG"), + string("NTG"), string("RHO") }; const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); @@ -561,9 +567,17 @@ void EclipseGridParser::convertToSI() key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" || key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" || key == "SGAS" || key == "SWAT" || key == "SOIL" || - key == "NTG" || key == "SWCR" || key == "SWL" || - key == "SWU" || key == "SOWCR" || key == "KRW" || - key == "KRWR" || key == "KRORW" || key == "KRO" || + key == "SWL" || key == "SWCR" || key == "SWU" || + key == "SGL " || key == "SGCR" || key == "SGU" || + key == "SOWCR" || key == "SOGCR" || key == "KRW" || + key == "KRWR" || key == "KRG" || key == "KRGR" || + key == "KRO" || key == "KRORW" || key == "KRORG" || + key == "ISWL" || key == "ISWCR" || key == "ISWU" || + key == "ISGL " || key == "ISGCR" || key == "ISGU" || + key == "ISOWCR" || key == "ISOGCR" || key == "IKRW" || + key == "IKRWR" || key == "IKRG" || key == "IKRGR" || + key == "IKRO" || key == "IKRORW" || key == "IKRORG" || + key == "NTG" || key == "RHO") /* nonstandard field with no unit logic. use with caution */ { unit = 1.0; do_convert = false; // Dimensionless keywords... @@ -571,8 +585,6 @@ void EclipseGridParser::convertToSI() unit = units_.pressure; } else if (key == "RS") { unit = units_.gasvol_s / units_.liqvol_s; - } else if (key == "RV") { - unit = units_.liqvol_s / units_.gasvol_s; } else if (key == "MAPAXES") { OPM_MESSAGE("Not applying units to MAPAXES yet!"); unit = 1.0; diff --git a/opm/core/io/eclipse/SpecialEclipseFields.hpp b/opm/core/io/eclipse/SpecialEclipseFields.hpp index 18b5635e..98432f4f 100644 --- a/opm/core/io/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/io/eclipse/SpecialEclipseFields.hpp @@ -2313,8 +2313,8 @@ struct ENPTVD : public SpecialBase { } virtual void read(std::istream & is) { - table_.resize(5); - std::vector > sub_table(5); + table_.resize(9); + std::vector > sub_table(9); while (!is.eof()) { if(is.peek() == int('/')) { if (sub_table[0].empty() && !(table_[0].empty())) { @@ -2332,19 +2332,14 @@ struct ENPTVD : public SpecialBase { if (data[0] == -1.0) { OPM_THROW(std::runtime_error, "Error reading ENPTVD data - depth can not be defaulted."); } - if ((data[4] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0) || (data[8] != -1.0)) { - OPM_THROW(std::runtime_error, "Error reading ENPTVD data - non-default values in column 5-7,9 not supported."); + for(std::vector >::size_type i=0; i> ignoreWhitespace; if(is.peek() == int('/')) { is >> ignoreLine; if (sub_table[0].size() >= 2) { - insertDefaultValues(sub_table, 5, -1.0, false); + insertDefaultValues(sub_table, 9, -1.0, false); std::vector >::iterator it_sub = sub_table.begin(); for(std::vector > >::size_type i=0; i >::size_type j=0; j::size_type k=0; k > >::size_type i=0; i > sub_table(5); + table_.resize(8); + std::vector > sub_table(8); while (!is.eof()) { if(is.peek() == int('/')) { if (sub_table[0].empty() && !(table_[0].empty())) { @@ -2419,19 +2414,15 @@ struct ENKRVD : public SpecialBase { if (data[0] == -1.0) { OPM_THROW(std::runtime_error, "Error reading ENKRVD data - depth can not be defaulted."); } - if ((data[2] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0)) { - OPM_THROW(std::runtime_error, "Error reading ENKRVD data - non-default values in column 3,6-7 not supported."); + + for(std::vector >::size_type i=0; i> ignoreWhitespace; if(is.peek() == int('/')) { is >> ignoreLine; if (sub_table[0].size() >= 2) { - insertDefaultValues(sub_table, 5, -1.0, false); + insertDefaultValues(sub_table, 8, -1.0, false); std::vector >::iterator it_sub = sub_table.begin(); for(std::vector > >::size_type i=0; i >::size_type j=0; j::size_type k=0; k > >::size_type i=0; i. +*/ +#ifndef SATFUNCBASE_HPP +#define SATFUNCBASE_HPP + +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + // Transforms for saturation table scaling + struct EPSTransforms { + struct Transform { + bool doNotScale; + bool do_3pt; + double smin; + double scr; + double sr; + double smax; + double slope1; + double slope2; + double scaleSat(double ss, double s_r, double s_cr, double s_max) const; + double scaleSatInv(double s, double s_r, double s_cr, double s_max) const; + double scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const; // Returns scaleSat'(s) + double scaleSatPc(double s, double s_min, double s_max) const; + double scaleSatDerivPc(double s, double s_min, double s_max) const; // Returns scaleSatPc'(s) + bool doKrMax; + bool doKrCrit; + bool doSatInterp; + double krsr; + double krmax; + double krSlopeMax; + 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)) + void printMe(std::ostream & out); + }; + + Transform wat; + Transform watoil; + Transform gas; + Transform gasoil; + }; + + // Hysteresis + struct SatHyst { + double sg_hyst; + double sg_shift; + double sow_hyst; + double sow_shift; + void printMe(std::ostream & out); + }; + + + template + class SatFuncBase : public BlackoilPhases + { + public: + void init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg, + const int samples); + void updateSatHyst(const double* s, + const EPSTransforms* epst, + const EPSTransforms* epst_hyst, + SatHyst* sat_hyst) const; + double smin_[PhaseUsage::MaxNumPhases]; + double smax_[PhaseUsage::MaxNumPhases]; + double krwmax_; // Max water relperm + double krgmax_; // Max gas relperm + double kromax_; // Max oil relperm + double swcr_; // Critical water saturation. + double sgcr_; // Critical gas saturation. + double krwr_; // Water relperm at critical oil-in-water saturation. + double krgr_; // Gas relperm at critical oil-in-gas saturation. + double sowcr_; // Critical oil-in-water saturation. + 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. + + protected: + PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. + TableType krw_; + TableType krow_; + TableType pcow_; + TableType krg_; + TableType krog_; + TableType pcog_; + double krocw_; // = krow_(s_wc) + + private: + void extendTable(const std::vector& xv, + std::vector& xv_ex, + double pm) const; + void initializeTableType(TableType& table, + const std::vector& arg, + const std::vector& value, + const int samples); + }; + + template + void SatFuncBase::init(const EclipseGridParser& deck, + const int table_num, + const PhaseUsage phase_usg, + const int samples) + { + phase_usage = phase_usg; + double swco = 0.0; + double swmax = 1.0; + if (phase_usage.phase_used[Aqua]) { + const SWOF::table_t& swof_table = deck.getSWOF().swof_; + const std::vector& sw = swof_table[table_num][0]; + const std::vector& krw = swof_table[table_num][1]; + const std::vector& krow = swof_table[table_num][2]; + const std::vector& pcow = swof_table[table_num][3]; + if (krw.front() != 0.0 || krow.back() != 0.0) { + OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)"); + } + + // Extend the tables with constant values such that the + // derivatives at the endpoints are zero + int n = sw.size(); + std::vector sw_ex(n+2); + std::vector krw_ex(n+2); + std::vector krow_ex(n+2); + std::vector pcow_ex(n+2); + + extendTable(sw,sw_ex,1); + extendTable(krw,krw_ex,0); + extendTable(krow,krow_ex,0); + extendTable(pcow,pcow_ex,0); + + initializeTableType(krw_,sw_ex, krw_ex, samples); + initializeTableType(krow_,sw_ex, krow_ex, samples); + initializeTableType(pcow_,sw_ex, pcow_ex, samples); + + krocw_ = krow[0]; // At connate water -> ecl. SWOF + swco = sw[0]; + smin_[phase_usage.phase_pos[Aqua]] = sw[0]; + swmax = sw.back(); + smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + + krwmax_ = krw.back(); + kromax_ = krow.front(); + swcr_ = swmax; + sowcr_ = 1.0 - swco; + krwr_ = krw.back(); + krorw_ = krow.front(); + for (std::vector::size_type i=1; i 0.0) { + swcr_ = sw[i-1]; + krorw_ = krow[i-1]; + break; + } + } + for (std::vector::size_type i=sw.size()-1; i>=1; --i) { + if (krow[i-1]> 0.0) { + sowcr_ = 1.0 - sw[i]; + krwr_ = krw[i]; + break; + } + } + } + if (phase_usage.phase_used[Vapour]) { + const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; + const std::vector& sg = sgof_table[table_num][0]; + const std::vector& krg = sgof_table[table_num][1]; + const std::vector& krog = sgof_table[table_num][2]; + const std::vector& pcog = sgof_table[table_num][3]; + + // Extend the tables with constant values such that the + // derivatives at the endpoints are zero + int n = sg.size(); + std::vector sg_ex(n+2); + std::vector krg_ex(n+2); + std::vector krog_ex(n+2); + std::vector pcog_ex(n+2); + + extendTable(sg,sg_ex,1); + extendTable(krg,krg_ex,0); + extendTable(krog,krog_ex,0); + extendTable(pcog,pcog_ex,0); + + initializeTableType(krg_,sg_ex, krg_ex, samples); + initializeTableType(krog_,sg_ex, krog_ex, samples); + initializeTableType(pcog_,sg_ex, pcog_ex, samples); + + smin_[phase_usage.phase_pos[Vapour]] = sg[0]; + if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { + OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << + ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); + } + smax_[phase_usage.phase_pos[Vapour]] = sg.back(); + smin_[phase_usage.phase_pos[Vapour]] = sg.front(); + krgmax_ = krg.back(); + + sgcr_ = sg.front(); + sogcr_ = 1.0 - sg.back(); + krgr_ = krg.back(); + krorg_ = krg.front(); + for (std::vector::size_type i=1; i 0.0) { + sgcr_ = sg[i-1]; + krorg_ = krog[i-1]; + break; + } + } + for (std::vector::size_type i=sg.size()-1; i>=1; --i) { + if (krog[i-1]> 0.0) { + sogcr_ = 1.0 - sg[i]; + krgr_ = krg[i]; + break; + } + } + + } + + if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) { + sowcr_ -= smin_[phase_usage.phase_pos[Vapour]]; + sogcr_ -= smin_[phase_usage.phase_pos[Aqua]]; + smin_[phase_usage.phase_pos[Liquid]] = 0.0; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]] + - smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ... + } else if (phase_usage.phase_used[Aqua]) { + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]]; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]; + } else if (phase_usage.phase_used[Vapour]) { + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]]; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]]; + } + } + + template + void SatFuncBase::updateSatHyst(const double* s, + const EPSTransforms* epst, + const EPSTransforms* epst_hyst, + SatHyst* sat_hyst) const + { + if (phase_usage.phase_used[Aqua] && phase_usage.phase_used[Vapour]) { //Water/Oil/Gas + int opos = phase_usage.phase_pos[BlackoilPhases::Liquid]; + int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour]; + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + if (s[opos] > sat_hyst->sow_hyst) + { + sat_hyst->sow_hyst = s[opos]; + double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]); + double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]); + sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst; + } + if (s[gpos] > sat_hyst->sg_hyst) + { + sat_hyst->sg_hyst = s[gpos]; + double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]); + double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]); + sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst; + } + } else if (phase_usage.phase_used[Aqua]) { //Water/oil + int opos = phase_usage.phase_pos[BlackoilPhases::Liquid]; + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + if (s[opos] > sat_hyst->sow_hyst) + { + sat_hyst->sow_hyst = s[opos]; + double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]); + double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]); + sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst; + } + } else if (phase_usage.phase_used[Vapour]) {//Gas/Oil + int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour]; + if (s[gpos] > sat_hyst->sg_hyst) + { + sat_hyst->sg_hyst = s[gpos]; + double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]); + double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]); + sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst; + } + } + } + + template + void SatFuncBase::extendTable(const std::vector& xv, + std::vector& xv_ex, + double pm) const + { + int n = xv.size(); + xv_ex[0] = xv[0]-pm; + xv_ex[n+1] = xv[n-1]+pm; + for (int i=0; i& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& pcow = swof_table[table_num][3]; - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sw.size(); - std::vector sw_ex(n+2); - std::vector krw_ex(n+2); - std::vector krow_ex(n+2); - std::vector pcow_ex(n+2); - - SatFuncGwsegUniform::ExtendTable(sw,sw_ex,1); - SatFuncGwsegUniform::ExtendTable(krw,krw_ex,0); - SatFuncGwsegUniform::ExtendTable(krow,krow_ex,0); - SatFuncGwsegUniform::ExtendTable(pcow,pcow_ex,0); - - buildUniformMonotoneTable(sw_ex, krw_ex, samples, krw_); - buildUniformMonotoneTable(sw_ex, krow_ex, samples, krow_); - buildUniformMonotoneTable(sw_ex, pcow_ex, samples, pcow_); - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage.phase_pos[Aqua]] = sw.back(); - } - if (phase_usage.phase_used[Vapour]) { - const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - const std::vector& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& pcog = sgof_table[table_num][3]; - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sg.size(); - std::vector sg_ex(n+2); - std::vector krg_ex(n+2); - std::vector krog_ex(n+2); - std::vector pcog_ex(n+2); - - SatFuncGwsegUniform::ExtendTable(sg,sg_ex,1); - SatFuncGwsegUniform::ExtendTable(krg,krg_ex,0); - SatFuncGwsegUniform::ExtendTable(krog,krog_ex,0); - SatFuncGwsegUniform::ExtendTable(pcog,pcog_ex,0); - - buildUniformMonotoneTable(sg_ex, krg_ex, samples, krg_); - buildUniformMonotoneTable(sg_ex, krog_ex, samples, krog_); - buildUniformMonotoneTable(sg_ex, pcog_ex, samples, pcog_); - smin_[phase_usage.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { - OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage.phase_pos[Vapour]] = sg.back(); - } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; - } - - void SatFuncGwsegUniform::ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const - { - int n = xv.size(); - xv_ex[0] = xv[0]-pm; - xv_ex[n+1] = xv[n-1]+pm; - for (int i=0; i& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& pcow = swof_table[table_num][3]; - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sw.size(); - std::vector sw_ex(n+2); - std::vector krw_ex(n+2); - std::vector krow_ex(n+2); - std::vector pcow_ex(n+2); - - SatFuncGwsegNonuniform::ExtendTable(sw,sw_ex,1); - SatFuncGwsegNonuniform::ExtendTable(krw,krw_ex,0); - SatFuncGwsegNonuniform::ExtendTable(krow,krow_ex,0); - SatFuncGwsegNonuniform::ExtendTable(pcow,pcow_ex,0); - - krw_ = NonuniformTableLinear(sw_ex, krw_ex); - krow_ = NonuniformTableLinear(sw_ex, krow_ex); - pcow_ = NonuniformTableLinear(sw_ex, pcow_ex); - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage.phase_pos[Aqua]] = sw.back(); - } - if (phase_usage.phase_used[Vapour]) { - const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - const std::vector& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& pcog = sgof_table[table_num][3]; - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sg.size(); - std::vector sg_ex(n+2); - std::vector krg_ex(n+2); - std::vector krog_ex(n+2); - std::vector pcog_ex(n+2); - - SatFuncGwsegNonuniform::ExtendTable(sg,sg_ex,1); - SatFuncGwsegNonuniform::ExtendTable(krg,krg_ex,0); - SatFuncGwsegNonuniform::ExtendTable(krog,krog_ex,0); - SatFuncGwsegNonuniform::ExtendTable(pcog,pcog_ex,0); - - - krg_ = NonuniformTableLinear(sg_ex, krg_ex); - krog_ = NonuniformTableLinear(sg_ex, krog_ex); - pcog_ = NonuniformTableLinear(sg_ex, pcog_ex); - smin_[phase_usage.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { - OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage.phase_pos[Vapour]] = sg.back(); - } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; - } - - void SatFuncGwsegNonuniform::ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const - { - int n = xv.size(); - xv_ex[0] = xv[0]-pm; - xv_ex[n+1] = xv[n-1]+pm; - for (int i=0; i -#include -#include -#include -#include +#include namespace Opm { - class SatFuncGwsegUniform : public BlackoilPhases + template + class SatFuncGwseg : public SatFuncBase { public: - void init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples); void evalKr(const double* s, double* kr) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalPc(const double* s, double* pc) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - void ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - double krwmax_; // Max water relperm - double kromax_; // Max oil relperm - double swcr_; // Critical water saturation. - double krwr_; // Water relperm at critical oil-in-water saturation. - double sowcr_; // Critical oil-in-water saturation. - double krorw_; // Oil relperm at critical water saturation. + + void evalKr(const double* s, double* kr, const EPSTransforms* epst) const; + void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const; + void evalPc(const double* s, double* pc, const EPSTransforms* epst) const; + void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const; + private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - UniformTableLinear krw_; - UniformTableLinear krow_; - UniformTableLinear pcow_; - UniformTableLinear krg_; - UniformTableLinear krog_; - UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) + }; + + typedef SatFuncGwseg > SatFuncGwsegUniform; + typedef SatFuncGwseg > SatFuncGwsegNonuniform; - - class SatFuncGwsegNonuniform : public BlackoilPhases + template + void SatFuncGwseg::evalKr(const double* s, double* kr) const { - public: - void init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples); - void evalKr(const double* s, double* kr) const; - void evalKrDeriv(const double* s, double* kr, double* dkrds) const; - void evalPc(const double* s, double* pc) const; - void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - void ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const; + if (this->phase_usage.num_phases == 3) { + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; + const double sw = std::max(s[BlackoilPhases::Aqua], swco); + const double sg = s[BlackoilPhases::Vapour]; + const double eps = 1e-5; + swco = std::min(swco,sw-eps); + + // xw and xg are the fractions occupied by water and gas zones. + const double ssg = sw - swco + sg; + const double xw = (sw - swco) / ssg; + const double xg = 1 - xw; + const double ssw = sg + sw; + + const double krw = this->krw_(sw); + const double krg = this->krg_(sg); + const double krow = this->krow_(ssw); + const double krog = this->krog_(ssg); + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = xw*krow + xg*krog; + return; + } + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double krow = this->krow_(sw); + kr[wpos] = krw; + kr[opos] = krow; + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double krog = this->krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + + template + void SatFuncGwseg::evalKr(const double* s, double* kr, const EPSTransforms* epst) const + { + if (this->phase_usage.num_phases == 3) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + + // TODO Also consider connate gas ... + double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; + double swco = epst->wat.smin; + const double sw = std::max(s[BlackoilPhases::Aqua], swco); + const double sg = s[BlackoilPhases::Vapour]; + const double eps = 1e-6; + swco = std::min(swco,sw-eps); + const double ssw = sg + sw; + const double ssg = std::max(sg + sw - swco, eps); + const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs). + + double ssow = 1.0-ssw; + double ssog = 1.0-ssg-swco; + double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + + const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_); + const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_); + const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_); + const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_); + + // xw and xg are the fractions occupied by water and gas zones. + const double xw = (sw - swco) / d; + const double xg = 1 - xw; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = xw*krow + xg*krog; + return; + } + OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ..."); + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double krow = this->krow_(sw); + kr[wpos] = krw; + kr[opos] = krow; + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double krog = this->krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + template + void SatFuncGwseg::evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const + { + if (this->phase_usage.num_phases == 3) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + + // TODO Consider connate gas ... + double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; + double swco = epst->wat.smin; + const double sw = std::max(s[BlackoilPhases::Aqua], swco); + const double sg = s[BlackoilPhases::Vapour]; + const double eps = 1e-6; + swco = std::min(swco,sw-eps); + const double ssw = sg + sw; + const double ssg = std::max(sg + sw - swco, eps); + const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs). + + double ssow = 1.0-ssw; + double ssog = 1.0-ssg-swco; + + // The code below corresponds to EHYSTR * 0 * * KR/ + // - wettability properties water>oil>gas. + // - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase. + // The imb-curve currently only differs from drainage curves via endpoint scaling ... + + // Water - use drainage curve only + double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_); + + // Gas + double krg; + if (sg >= sat_hyst->sg_hyst) { // Drainage + double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_); + } else { // Imbibition + double sg_shifted = sg + sat_hyst->sg_shift; + double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_); + } + + // Oil in water + double krow; + if (ssow >= sat_hyst->sow_hyst) { // Drainage + double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_); + } else { // Imbibition + double ssow_shifted = ssow + sat_hyst->sow_shift; + double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_); + } + + // Oil in gas and connate water - use drainage curve only + double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_); + + // xw and xg are the fractions occupied by water and gas zones. + const double xw = (sw - swco) / d; + const double xg = 1 - xw; + + // relperms + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = xw*krow + xg*krog; + + return; + } + OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ..."); + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double krow = this->krow_(sw); + kr[wpos] = krw; + kr[opos] = krow; + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double krog = this->krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + template + void SatFuncGwseg::evalKrDeriv(const double* s, double* kr, double* dkrds) const + { + const int np = this->phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + + double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; + const double sw = std::max(s[BlackoilPhases::Aqua], swco); + const double sg = s[BlackoilPhases::Vapour]; + const double eps = 1e-5; + swco = std::min(swco,sw-eps); + + // xw and xg are the fractions occupied by water and gas zones. + const double ssw = sg + sw; + // d = ssg = sw - swco + sg (using 'd' for consistency with mrst docs). + const double d = sg + sw - swco; + const double xw = (sw - swco) / d; + + const double krw = this->krw_(sw); + const double krg = this->krg_(sg); + const double krow = this->krow_(ssw); + const double krog = this->krog_(d); + + const double xg = 1 - xw; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = xw*krow + xg*krog; + + // Derivatives. + const double dkrww = this->krw_.derivative(sw); + const double dkrgg = this->krg_.derivative(sg); + const double dkrow = this->krow_.derivative(ssw); + const double dkrog = this->krog_.derivative(d); + dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog; + dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg; + return; + } + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double krow = this->krow_(sw); + double dkrow = this->krow_.derivative(sw); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + + } + + template + void SatFuncGwseg::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const + { + const int np = this->phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + + // TODO Also consider connate gas ... + double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; + double swco = epst->wat.smin; + const double sw = std::max(s[BlackoilPhases::Aqua], swco); + const double sg = s[BlackoilPhases::Vapour]; + const double eps = 1e-6; + swco = std::min(swco,sw-eps); + const double ssw = sg + sw; + const double ssg = std::max(sg + sw - swco, eps); + const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs). + + double ssow = 1.0-ssw; + double ssog = 1.0-ssg-swco; + double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + + const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_); + const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_); + const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_); + const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_); + + // xw and xg are the fractions occupied by water and gas zones. + const double xw = (sw - swco) / d; + const double xg = 1 - xw; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = xw*krow + xg*krog; + + // Derivatives. + double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw)); + double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg)); + double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow)); + double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco)); + dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog; + dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg; + return; + } + OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ..."); + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double krow = this->krow_(sw); + double dkrow = this->krow_.derivative(sw); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + } + + template + void SatFuncGwseg::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const + { + const int np = this->phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + + // Relative permeability model based on segregation of water + // and gas, with oil present in both water and gas zones. + + // TODO Also consider connate gas ... + double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]]; + double swco = epst->wat.smin; + const double sw = std::max(s[BlackoilPhases::Aqua], swco); + const double sg = s[BlackoilPhases::Vapour]; + const double eps = 1e-6; + swco = std::min(swco,sw-eps); + const double ssw = sg + sw; + const double ssg = std::max(sg + sw - swco, eps); + const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs). + + double ssow = 1.0-ssw; + double ssog = 1.0-ssg-swco; + + // The code below corresponds to EHYSTR * 0 * * KR/ + // - wettability properties water>oil>gas. + // - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase. + // The imb-curve currently only differs from drainage curves via endpoint scaling ... + + // Water - use drainage curve only + double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_); + double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw)); + + // Gas + double krg, dkrgg; + if (sg >= sat_hyst->sg_hyst) { // Drainage + double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_); + dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg)); + } else { // Imbibition + double sg_shifted = sg + sat_hyst->sg_shift; + double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _dsdsg = epst_hyst->gas.scaleSatDeriv(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_); + dkrgg = _dsdsg*epst_hyst->gas.scaleKrDeriv(sg_shifted, this->krg_.derivative(_sg)); + } + + // Oil in water + double krow, dkrow; + if (ssow >= sat_hyst->sow_hyst) { // Drainage + double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_); + dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow)); + } else { // Imbibition + double ssow_shifted = ssow + sat_hyst->sow_shift; + double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _dsdssow = epst_hyst->watoil.scaleSatDeriv(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_); + dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(1.0-_ssow)); + } + + // Oil in gas and connate water - use drainage curve only + double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_); + double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco)); + + // xw and xg are the fractions occupied by water and gas zones. + const double xw = (sw - swco) / d; + const double xg = 1 - xw; + + // relperms + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = xw*krow + xg*krog; + + // Derivatives. + dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog; + dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg; + return; + } + OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ..."); + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double krow = this->krow_(sw); + double dkrow = this->krow_.derivative(sw); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + } + + template + void SatFuncGwseg::evalPc(const double* s, double* pc) const + { + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + pc[pos] = this->pcow_(s[pos]); + } + if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + pc[pos] = this->pcog_(s[pos]); + } + } + + + template + void SatFuncGwseg::evalPc(const double* s, double* pc, const EPSTransforms* epst) const + { + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + 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); + } + 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); + } + } + + template + void SatFuncGwseg::evalPcDeriv(const double* s, double* pc, double* dpcds) const + { + // The problem of determining three-phase capillary pressures + // is very hard experimentally, usually one extends two-phase + // data (as for relative permeability). + // In our approach the derivative matrix is quite sparse, only + // the diagonal elements corresponding to non-oil phases are + // (potentially) nonzero. + const int np = this->phase_usage.num_phases; + std::fill(dpcds, dpcds + np*np, 0.0); + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + pc[pos] = this->pcow_(s[pos]); + dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]); + } + if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + pc[pos] = this->pcog_(s[pos]); + dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]); + } + } + + + template + void SatFuncGwseg::evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const + { + // The problem of determining three-phase capillary pressures + // is very hard experimentally, usually one extends two-phase + // data (as for relative permeability). + // In our approach the derivative matrix is quite sparse, only + // the diagonal elements corresponding to non-oil phases are + // (potentially) nonzero. + const int np = this->phase_usage.num_phases; + std::fill(dpcds, dpcds + np*np, 0.0); + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + 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]); + double _dsdsw = epst->wat.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]); + dpcds[np*pos + pos] = _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); + double _dsdsg = epst->gas.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]); + dpcds[np*pos + pos] = _dsdsg*this->pcog_.derivative(_sg); + } + } - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - double krwmax_; // Max water relperm - double kromax_; // Max oil relperm - double swcr_; // Critical water saturation. - double krwr_; // Water relperm at critical oil-in-water saturation. - double sowcr_; // Critical oil-in-water saturation. - double krorw_; // Oil relperm at critical water saturation. - private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - NonuniformTableLinear krw_; - NonuniformTableLinear krow_; - NonuniformTableLinear pcow_; - NonuniformTableLinear krg_; - NonuniformTableLinear krog_; - NonuniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - }; } // namespace Opm #endif // SATFUNCGWSEG_HPP diff --git a/opm/core/props/satfunc/SatFuncSimple.cpp b/opm/core/props/satfunc/SatFuncSimple.cpp index 58ea618d..940f12de 100644 --- a/opm/core/props/satfunc/SatFuncSimple.cpp +++ b/opm/core/props/satfunc/SatFuncSimple.cpp @@ -30,536 +30,194 @@ namespace Opm { - - - - void SatFuncSimpleUniform::init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples) + // SatFuncSimple.cpp can now be removed and the code below moved to SatFuncBase.cpp ... + + template<> + void SatFuncBase >::initializeTableType(NonuniformTableLinear & table, + const std::vector& arg, + const std::vector& value, + const int samples) { - phase_usage = phase_usg; - double swco = 0.0; - double swmax = 1.0; - if (phase_usage.phase_used[Aqua]) { - const SWOF::table_t& swof_table = deck.getSWOF().swof_; - const std::vector& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& pcow = swof_table[table_num][3]; - if (krw.front() != 0.0 || krow.back() != 0.0) { - OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)"); - } - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sw.size(); - std::vector sw_ex(n+2); - std::vector krw_ex(n+2); - std::vector krow_ex(n+2); - std::vector pcow_ex(n+2); - - SatFuncSimpleUniform::ExtendTable(sw,sw_ex,1); - SatFuncSimpleUniform::ExtendTable(krw,krw_ex,0); - SatFuncSimpleUniform::ExtendTable(krow,krow_ex,0); - SatFuncSimpleUniform::ExtendTable(pcow,pcow_ex,0); - - buildUniformMonotoneTable(sw_ex, krw_ex, samples, krw_); - buildUniformMonotoneTable(sw_ex, krow_ex, samples, krow_); - buildUniformMonotoneTable(sw_ex, pcow_ex, samples, pcow_); - - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage.phase_pos[Aqua]] = sw.back(); - - krwmax_ = krw.back(); - kromax_ = krow.front(); - swcr_ = swmax; - sowcr_ = 1.0 - swco; - krwr_ = krw.back(); - krorw_ = krow.front(); - for (std::vector::size_type i=1; i 0.0) { - swcr_ = sw[i-1]; - krorw_ = krow[i-1]; - break; - } - } - for (std::vector::size_type i=sw.size()-1; i>=1; --i) { - if (krow[i-1]> 0.0) { - sowcr_ = 1.0 - sw[i]; - krwr_ = krw[i]; - break; - } - } - } - if (phase_usage.phase_used[Vapour]) { - const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - const std::vector& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& pcog = sgof_table[table_num][3]; - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sg.size(); - std::vector sg_ex(n+2); - std::vector krg_ex(n+2); - std::vector krog_ex(n+2); - std::vector pcog_ex(n+2); - - SatFuncSimpleUniform::ExtendTable(sg,sg_ex,1); - SatFuncSimpleUniform::ExtendTable(krg,krg_ex,0); - SatFuncSimpleUniform::ExtendTable(krog,krog_ex,0); - SatFuncSimpleUniform::ExtendTable(pcog,pcog_ex,0); - - buildUniformMonotoneTable(sg_ex, krg_ex, samples, krg_); - buildUniformMonotoneTable(sg_ex, krog_ex, samples, krog_); - buildUniformMonotoneTable(sg_ex, pcog_ex, samples, pcog_); - smin_[phase_usage.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { - OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage.phase_pos[Vapour]] = sg.back(); - } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; + table = NonuniformTableLinear(arg, value); } - - - void SatFuncSimpleUniform::ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const + template<> + void SatFuncBase >::initializeTableType(UniformTableLinear & table, + const std::vector& arg, + const std::vector& value, + const int samples) { - int n = xv.size(); - xv_ex[0] = xv[0]-pm; - xv_ex[n+1] = xv[n-1]+pm; - for (int i=0; i= smax) ? s_max : s_cr + (s-scr)*slope1; } - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double so = s[opos]; - double krow = krow_(1.0-so); - kr[wpos] = krw; - kr[opos] = krow; + } else if (s <= sr) { + return (s <= scr) ? s_cr : s_cr+(s-scr)*slope1; } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double krog = krog_(sg); - kr[gpos] = krg; - kr[opos] = krog; + return (s >= smax) ? s_max : s_r+(s-sr)*slope2; } } - - void SatFuncSimpleUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const + double EPSTransforms::Transform::scaleSatInv(double ss, double s_r, double s_cr, double s_max) const { - const int np = phase_usage.num_phases; - std::fill(dkrds, dkrds + np*np, 0.0); - - if (np == 3) { - // A simplified relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krow = krow_(sw + sg); - double dkrow = krow_.derivative(sw + sg); - // double krog = krog_(sg); - // double dkrog = krog_.derivative(sg); - // double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krow; - //krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; + if (doNotScale) { + return ss; + } else if (!do_3pt) { // 2-pt + if (ss <= s_cr) { + return scr; + } else { + return (ss >= s_max) ? smax : scr + (ss-s_cr)/slope1; } - dkrds[Aqua + Aqua*np] = dkrww; - dkrds[Vapour + Vapour*np] = dkrgg; - //dkrds[Liquid + Aqua*np] = dkrow; - dkrds[Liquid + Liquid*np] = -dkrow; - //krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); - dkrds[Liquid + Vapour*np] = 0.0; - //krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) - //+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double so = s[opos]; - double krow = krow_(1.0-so); - double dkrow = krow_.derivative(1.0-so); - kr[wpos] = krw; - kr[opos] = krow; - dkrds[wpos + wpos*np] = dkrww; - dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else if (ss <= s_r) { + return (ss <= s_cr) ? scr : scr+(ss-s_cr)/slope1; } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krog = krog_(sg); - double dkrog = krog_.derivative(sg); - kr[gpos] = krg; - kr[opos] = krog; - dkrds[gpos + gpos*np] = dkrgg; - dkrds[opos + gpos*np] = dkrog; - } - - } - - - void SatFuncSimpleUniform::evalPc(const double* s, double* pc) const - { - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); + return (ss >= s_max) ? smax : sr+(ss-s_r)/slope2; } } - void SatFuncSimpleUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const + double EPSTransforms::Transform::scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const { - // The problem of determining three-phase capillary pressures - // is very hard experimentally, usually one extends two-phase - // data (as for relative permeability). - // In our approach the derivative matrix is quite sparse, only - // the diagonal elements corresponding to non-oil phases are - // (potentially) nonzero. - const int np = phase_usage.num_phases; - std::fill(dpcds, dpcds + np*np, 0.0); - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - dpcds[np*pos + pos] = pcow_.derivative(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); - dpcds[np*pos + pos] = pcog_.derivative(s[pos]); - } - } - - - - - - // ====== Methods for SatFuncSimpleNonuniform ====== - - - - void SatFuncSimpleNonuniform::ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const - { - int n = xv.size(); - xv_ex[0] = xv[0]-pm; - xv_ex[n+1] = xv[n-1]+pm; - for (int i=0; i& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& pcow = swof_table[table_num][3]; - if (krw.front() != 0.0 || krow.back() != 0.0) { - OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)"); + if (doNotScale) { + return 1.0; + } else if (!do_3pt) { // 2-pt + if (s <= scr) { + return 0.0; + } else { + return (s >= smax) ? 0.0 : slope1; } - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sw.size(); - std::vector sw_ex(n+2); - std::vector krw_ex(n+2); - std::vector krow_ex(n+2); - std::vector pcow_ex(n+2); - - SatFuncSimpleNonuniform::ExtendTable(sw,sw_ex,1); - SatFuncSimpleNonuniform::ExtendTable(krw,krw_ex,0); - SatFuncSimpleNonuniform::ExtendTable(krow,krow_ex,0); - SatFuncSimpleNonuniform::ExtendTable(pcow,pcow_ex,0); - - krw_ = NonuniformTableLinear(sw_ex, krw_ex); - krow_ = NonuniformTableLinear(sw_ex, krow_ex); - pcow_ = NonuniformTableLinear(sw_ex, pcow_ex); - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage.phase_pos[Aqua]] = sw.back(); - - krwmax_ = krw.back(); - kromax_ = krow.front(); - swcr_ = swmax; - sowcr_ = 1.0 - swco; - krwr_ = krw.back(); - krorw_ = krow.front(); - for (std::vector::size_type i=1; i 0.0) { - swcr_ = sw[i-1]; - krorw_ = krow[i-1]; - break; - } - } - for (std::vector::size_type i=sw.size()-1; i>=1; --i) { - if (krow[i-1]> 0.0) { - sowcr_ = 1.0 - sw[i]; - krwr_ = krw[i]; - break; - } - } - - } - if (phase_usage.phase_used[Vapour]) { - const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - const std::vector& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& pcog = sgof_table[table_num][3]; - - // Extend the tables with constant values such that the - // derivatives at the endpoints are zero - int n = sg.size(); - std::vector sg_ex(n+2); - std::vector krg_ex(n+2); - std::vector krog_ex(n+2); - std::vector pcog_ex(n+2); - - SatFuncSimpleNonuniform::ExtendTable(sg,sg_ex,1); - SatFuncSimpleNonuniform::ExtendTable(krg,krg_ex,0); - SatFuncSimpleNonuniform::ExtendTable(krog,krog_ex,0); - SatFuncSimpleNonuniform::ExtendTable(pcog,pcog_ex,0); - - krg_ = NonuniformTableLinear(sg_ex, krg_ex); - krog_ = NonuniformTableLinear(sg_ex, krog_ex); - pcog_ = NonuniformTableLinear(sg_ex, pcog_ex); - smin_[phase_usage.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { - OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage.phase_pos[Vapour]] = sg.back(); - } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; - } - - - void SatFuncSimpleNonuniform::evalKr(const double* s, double* kr) const - { - if (phase_usage.num_phases == 3) { - // A simplified relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double krg = krg_(sg); - double krow = krow_(sw + sg); // = 1 - so - // double krog = krog_(sg); // = 1 - so - sw - // double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krow; - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; - } - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double so = s[opos]; - double krow = krow_(1.0-so); - kr[wpos] = krw; - kr[opos] = krow; + } else if (s <= sr) { + return (s <= scr) ? 0.0 : slope1; } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double krog = krog_(sg); - kr[gpos] = krg; - kr[opos] = krog; + return (s >= smax) ? 0.0 : slope2; } } - void SatFuncSimpleNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const + double EPSTransforms::Transform::scaleSatPc(double s, double s_min, double s_max) const { - const int np = phase_usage.num_phases; - std::fill(dkrds, dkrds + np*np, 0.0); - - if (np == 3) { - // A simplified relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krow = krow_(sw + sg); - double dkrow = krow_.derivative(sw + sg); - // double krog = krog_(sg); - // double dkrog = krog_.derivative(sg); - // double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krow; - //krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; - } - dkrds[Aqua + Aqua*np] = dkrww; - dkrds[Vapour + Vapour*np] = dkrgg; - //dkrds[Liquid + Aqua*np] = dkrow; - dkrds[Liquid + Liquid*np] = -dkrow; - //krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); - dkrds[Liquid + Vapour*np] = 0.0; - //krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) - //+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double so = s[opos]; - double krow = krow_(1.0-so); - double dkrow = krow_.derivative(1.0-so); - kr[wpos] = krw; - kr[opos] = krow; - dkrds[wpos + wpos*np] = dkrww; - dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + if (doNotScale) { + return s; + } else if (s<=smin) { + return s_min; + } else if (s <= smax) { + return s_min + (s-smin)*(s_max-s_min)/(smax-smin); } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krog = krog_(sg); - double dkrog = krog_.derivative(sg); - kr[gpos] = krg; - kr[opos] = krog; - dkrds[gpos + gpos*np] = dkrgg; - dkrds[opos + gpos*np] = dkrog; + return s_max; } - } - - - void SatFuncSimpleNonuniform::evalPc(const double* s, double* pc) const + double EPSTransforms::Transform::scaleSatDerivPc(double s, double s_min, double s_max) const { - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); + if (doNotScale) { + return 1.0; + } else if (s -#include -#include -#include -#include +#include namespace Opm { - class SatFuncSimpleUniform : public BlackoilPhases + template + class SatFuncSimple : public SatFuncBase { public: - void init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples); void evalKr(const double* s, double* kr) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalPc(const double* s, double* pc) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - void ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - double krwmax_; // Max water relperm - double kromax_; // Max oil relperm - double swcr_; // Critical water saturation. - double krwr_; // Water relperm at critical oil-in-water saturation. - double sowcr_; // Critical oil-in-water saturation. - double krorw_; // Oil relperm at critical water saturation. + + void evalKr(const double* s, double* kr, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} + void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const + {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const; + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const + {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} + void evalPc(const double* s, double* pc, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} + void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} + private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - UniformTableLinear krw_; - UniformTableLinear krow_; - UniformTableLinear pcow_; - UniformTableLinear krg_; - UniformTableLinear krog_; - UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) + }; + + typedef SatFuncSimple > SatFuncSimpleUniform; + typedef SatFuncSimple > SatFuncSimpleNonuniform; - - class SatFuncSimpleNonuniform : public BlackoilPhases + template + void SatFuncSimple::evalKr(const double* s, double* kr) const { - public: - void init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples); - void evalKr(const double* s, double* kr) const; - void evalKrDeriv(const double* s, double* kr, double* dkrds) const; - void evalPc(const double* s, double* pc) const; - void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - void ExtendTable(const std::vector& xv, - std::vector& xv_ex, - double pm) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - double krwmax_; // Max water relperm - double kromax_; // Max oil relperm - double swcr_; // Critical water saturation. - double krwr_; // Water relperm at critical oil-in-water saturation. - double sowcr_; // Critical oil-in-water saturation. - double krorw_; // Oil relperm at critical water saturation. - private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - NonuniformTableLinear krw_; - NonuniformTableLinear krow_; - NonuniformTableLinear pcow_; - NonuniformTableLinear krg_; - NonuniformTableLinear krog_; - NonuniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - }; + if (this->phase_usage.num_phases == 3) { + // A simplified relative permeability model. + double sw = s[BlackoilPhases::Aqua]; + double sg = s[BlackoilPhases::Vapour]; + double krw = this->krw_(sw); + double krg = this->krg_(sg); + double krow = this->krow_(sw + sg); // = 1 - so + // double krog = krog_(sg); // = 1 - so - sw + // double krocw = krocw_; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = krow; + if (kr[BlackoilPhases::Liquid] < 0.0) { + kr[BlackoilPhases::Liquid] = 0.0; + } + return; + } + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double so = s[opos]; + double krow = this->krow_(1.0-so); + kr[wpos] = krw; + kr[opos] = krow; + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double krog = this->krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + template + void SatFuncSimple::evalKrDeriv(const double* s, double* kr, double* dkrds) const + { + const int np = this->phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + // A simplified relative permeability model. + double sw = s[BlackoilPhases::Aqua]; + double sg = s[BlackoilPhases::Vapour]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krow = this->krow_(sw + sg); + double dkrow = this->krow_.derivative(sw + sg); + // double krog = krog_(sg); + // double dkrog = krog_.derivative(sg); + // double krocw = krocw_; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = krow; + //krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); + if (kr[BlackoilPhases::Liquid] < 0.0) { + kr[BlackoilPhases::Liquid] = 0.0; + } + dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww; + dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg; + //dkrds[Liquid + Aqua*np] = dkrow; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow; + //krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = 0.0; + //krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) + //+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); + return; + } + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double so = s[opos]; + double krow = this->krow_(1.0-so); + double dkrow = this->krow_.derivative(1.0-so); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + + } + + template + void SatFuncSimple::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const + { + const int np = this->phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + // A simplified relative permeability model. + // Define KR(s) = scaleKr(kr(scalSat(s))) + // Thus KR'(s) = scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))*scaleSat'(s) + double _sw = epst->wat.scaleSat(s[wpos], 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double _dsdsw = epst->wat.scaleSatDeriv(s[wpos], 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]); + double _sg = epst->gas.scaleSat(s[gpos], 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _dsdsg = epst->gas.scaleSatDeriv(s[gpos], 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]); + double _sow = epst->watoil.scaleSat(1.0-s[wpos]-s[gpos], 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + double _dsdsow = epst->watoil.scaleSatDeriv(1.0-s[wpos]-s[gpos], 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + //double _sog = epst->gasoil.scaleSat(1.0-s[wpos]-s[gpos], 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + //double _dsdsog = epst->gasoil.scaleSatDeriv(1.0-s[wpos]-s[gpos], 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]); + + double krw = epst->wat.scaleKr(s[wpos], this->krw_(_sw), this->krwr_); + double dkrww = _dsdsw*epst->wat.scaleKrDeriv(s[wpos], this->krw_.derivative(_sw)); + double krg = epst->gas.scaleKr(s[gpos], this->krg_(_sg), this->krgr_); + double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(s[gpos], this->krg_.derivative(_sg)); + // TODO Check the arguments to the krow- and krog-tables below... + double krow = epst->watoil.scaleKr(1.0-s[wpos]-s[gpos], this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ???? + double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ???? + //double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-s[wpos]-s[gpos], this->krorg_); // ???? + //double dkrog = _dsdsog*epst->gasoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krog_.derivative(1.0-_sog-this->smin_[wpos])); // ???? + // double krocw = krocw_; + kr[wpos] = krw; + kr[gpos] = krg; + kr[BlackoilPhases::Liquid] = krow; + //krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); + if (kr[BlackoilPhases::Liquid] < 0.0) { + kr[BlackoilPhases::Liquid] = 0.0; + } + dkrds[wpos + wpos*np] = dkrww; + dkrds[gpos + gpos*np] = dkrgg; + //dkrds[Liquid + Aqua*np] = dkrow; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow; + //krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); + dkrds[BlackoilPhases::Liquid + gpos*np] = 0.0; + //krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) + //+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); + return; + } + OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ..."); + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double so = s[opos]; + double krow = this->krow_(1.0-so); + double dkrow = this->krow_.derivative(1.0-so); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + + } + + template + void SatFuncSimple::evalPc(const double* s, double* pc) const + { + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + pc[pos] = this->pcow_(s[pos]); + } + if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + pc[pos] = this->pcog_(s[pos]); + } + } + + template + void SatFuncSimple::evalPcDeriv(const double* s, double* pc, double* dpcds) const + { + // The problem of determining three-phase capillary pressures + // is very hard experimentally, usually one extends two-phase + // data (as for relative permeability). + // In our approach the derivative matrix is quite sparse, only + // the diagonal elements corresponding to non-oil phases are + // (potentially) nonzero. + const int np = this->phase_usage.num_phases; + std::fill(dpcds, dpcds + np*np, 0.0); + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + pc[pos] = this->pcow_(s[pos]); + dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]); + } + if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + pc[pos] = this->pcog_(s[pos]); + dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]); + } + } } // namespace Opm #endif // SATFUNCSIMPLE_HPP diff --git a/opm/core/props/satfunc/SatFuncStone2.cpp b/opm/core/props/satfunc/SatFuncStone2.cpp index 670616be..a9857a92 100644 --- a/opm/core/props/satfunc/SatFuncStone2.cpp +++ b/opm/core/props/satfunc/SatFuncStone2.cpp @@ -29,390 +29,5 @@ namespace Opm { - - - - - void SatFuncStone2Uniform::init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples) - { - phase_usage = phase_usg; - double swco = 0.0; - double swmax = 1.0; - if (phase_usage.phase_used[Aqua]) { - const SWOF::table_t& swof_table = deck.getSWOF().swof_; - const std::vector& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& pcow = swof_table[table_num][3]; - buildUniformMonotoneTable(sw, krw, samples, krw_); - buildUniformMonotoneTable(sw, krow, samples, krow_); - buildUniformMonotoneTable(sw, pcow, samples, pcow_); - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage.phase_pos[Aqua]] = sw.back(); - } - if (phase_usage.phase_used[Vapour]) { - const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - const std::vector& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& pcog = sgof_table[table_num][3]; - buildUniformMonotoneTable(sg, krg, samples, krg_); - buildUniformMonotoneTable(sg, krog, samples, krog_); - buildUniformMonotoneTable(sg, pcog, samples, pcog_); - smin_[phase_usage.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { - OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage.phase_pos[Vapour]] = sg.back(); - } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; - } - - - void SatFuncStone2Uniform::evalKr(const double* s, double* kr) const - { - if (phase_usage.num_phases == 3) { - // Stone-II relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double krg = krg_(sg); - double krow = krow_(sw + sg); // = 1 - so - double krog = krog_(sg); // = 1 - so - sw - double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; - } - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double krow = krow_(sw); - kr[wpos] = krw; - kr[opos] = krow; - } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double krog = krog_(sg); - kr[gpos] = krg; - kr[opos] = krog; - } - } - - - void SatFuncStone2Uniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const - { - const int np = phase_usage.num_phases; - std::fill(dkrds, dkrds + np*np, 0.0); - - if (np == 3) { - // Stone-II relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krow = krow_(sw + sg); - double dkrow = krow_.derivative(sw + sg); - double krog = krog_(sg); - double dkrog = krog_.derivative(sg); - double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; - } - dkrds[Aqua + Aqua*np] = dkrww; - dkrds[Vapour + Vapour*np] = dkrgg; - dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); - dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) - + krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double krow = krow_(sw); - double dkrow = krow_.derivative(sw); - kr[wpos] = krw; - kr[opos] = krow; - dkrds[wpos + wpos*np] = dkrww; - dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. - } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krog = krog_(sg); - double dkrog = krog_.derivative(sg); - kr[gpos] = krg; - kr[opos] = krog; - dkrds[gpos + gpos*np] = dkrgg; - dkrds[opos + gpos*np] = dkrog; - } - - } - - - void SatFuncStone2Uniform::evalPc(const double* s, double* pc) const - { - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); - } - } - - void SatFuncStone2Uniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const - { - // The problem of determining three-phase capillary pressures - // is very hard experimentally, usually one extends two-phase - // data (as for relative permeability). - // In our approach the derivative matrix is quite sparse, only - // the diagonal elements corresponding to non-oil phases are - // (potentially) nonzero. - const int np = phase_usage.num_phases; - std::fill(dpcds, dpcds + np*np, 0.0); - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - dpcds[np*pos + pos] = pcow_.derivative(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); - dpcds[np*pos + pos] = pcog_.derivative(s[pos]); - } - } - - - - - - - // ====== Methods for SatFuncStone2Nonuniform ====== - - - - - void SatFuncStone2Nonuniform::init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int /*samples*/) - { - phase_usage = phase_usg; - double swco = 0.0; - double swmax = 1.0; - if (phase_usage.phase_used[Aqua]) { - const SWOF::table_t& swof_table = deck.getSWOF().swof_; - const std::vector& sw = swof_table[table_num][0]; - const std::vector& krw = swof_table[table_num][1]; - const std::vector& krow = swof_table[table_num][2]; - const std::vector& pcow = swof_table[table_num][3]; - krw_ = NonuniformTableLinear(sw, krw); - krow_ = NonuniformTableLinear(sw, krow); - pcow_ = NonuniformTableLinear(sw, pcow); - krocw_ = krow[0]; // At connate water -> ecl. SWOF - swco = sw[0]; - smin_[phase_usage.phase_pos[Aqua]] = sw[0]; - swmax = sw.back(); - smax_[phase_usage.phase_pos[Aqua]] = sw.back(); - } - if (phase_usage.phase_used[Vapour]) { - const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; - const std::vector& sg = sgof_table[table_num][0]; - const std::vector& krg = sgof_table[table_num][1]; - const std::vector& krog = sgof_table[table_num][2]; - const std::vector& pcog = sgof_table[table_num][3]; - krg_ = NonuniformTableLinear(sg, krg); - krog_ = NonuniformTableLinear(sg, krog); - pcog_ = NonuniformTableLinear(sg, pcog); - smin_[phase_usage.phase_pos[Vapour]] = sg[0]; - if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { - OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << - ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); - } - smax_[phase_usage.phase_pos[Vapour]] = sg.back(); - } - // These only consider water min/max sats. Consider gas sats? - smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax; - smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco; - } - - - void SatFuncStone2Nonuniform::evalKr(const double* s, double* kr) const - { - if (phase_usage.num_phases == 3) { - // Stone-II relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double krg = krg_(sg); - double krow = krow_(sw + sg); // = 1 - so - double krog = krog_(sg); // = 1 - so - sw - double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; - } - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double krow = krow_(sw); - kr[wpos] = krw; - kr[opos] = krow; - } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double krog = krog_(sg); - kr[gpos] = krg; - kr[opos] = krog; - } - } - - - void SatFuncStone2Nonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const - { - const int np = phase_usage.num_phases; - std::fill(dkrds, dkrds + np*np, 0.0); - - if (np == 3) { - // Stone-II relative permeability model. - double sw = s[Aqua]; - double sg = s[Vapour]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krow = krow_(sw + sg); - double dkrow = krow_.derivative(sw + sg); - double krog = krog_(sg); - double dkrog = krog_.derivative(sg); - double krocw = krocw_; - kr[Aqua] = krw; - kr[Vapour] = krg; - kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); - if (kr[Liquid] < 0.0) { - kr[Liquid] = 0.0; - } - dkrds[Aqua + Aqua*np] = dkrww; - dkrds[Vapour + Vapour*np] = dkrgg; - dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); - dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) - + krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); - return; - } - // We have a two-phase situation. We know that oil is active. - if (phase_usage.phase_used[Aqua]) { - int wpos = phase_usage.phase_pos[Aqua]; - int opos = phase_usage.phase_pos[Liquid]; - double sw = s[wpos]; - double krw = krw_(sw); - double dkrww = krw_.derivative(sw); - double krow = krow_(sw); - double dkrow = krow_.derivative(sw); - kr[wpos] = krw; - kr[opos] = krow; - dkrds[wpos + wpos*np] = dkrww; - dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. - } else { - assert(phase_usage.phase_used[Vapour]); - int gpos = phase_usage.phase_pos[Vapour]; - int opos = phase_usage.phase_pos[Liquid]; - double sg = s[gpos]; - double krg = krg_(sg); - double dkrgg = krg_.derivative(sg); - double krog = krog_(sg); - double dkrog = krog_.derivative(sg); - kr[gpos] = krg; - kr[opos] = krog; - dkrds[gpos + gpos*np] = dkrgg; - dkrds[opos + gpos*np] = dkrog; - } - - } - - - void SatFuncStone2Nonuniform::evalPc(const double* s, double* pc) const - { - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); - } - } - - void SatFuncStone2Nonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const - { - // The problem of determining three-phase capillary pressures - // is very hard experimentally, usually one extends two-phase - // data (as for relative permeability). - // In our approach the derivative matrix is quite sparse, only - // the diagonal elements corresponding to non-oil phases are - // (potentially) nonzero. - const int np = phase_usage.num_phases; - std::fill(dpcds, dpcds + np*np, 0.0); - pc[phase_usage.phase_pos[Liquid]] = 0.0; - if (phase_usage.phase_used[Aqua]) { - int pos = phase_usage.phase_pos[Aqua]; - pc[pos] = pcow_(s[pos]); - dpcds[np*pos + pos] = pcow_.derivative(s[pos]); - } - if (phase_usage.phase_used[Vapour]) { - int pos = phase_usage.phase_pos[Vapour]; - pc[pos] = pcog_(s[pos]); - dpcds[np*pos + pos] = pcog_.derivative(s[pos]); - } - } - - - - - -} // namespace Opm +// SatFuncStone2.cpp can now be removed +} // namespace Opm diff --git a/opm/core/props/satfunc/SatFuncStone2.hpp b/opm/core/props/satfunc/SatFuncStone2.hpp index c698948a..4eb5e2fd 100644 --- a/opm/core/props/satfunc/SatFuncStone2.hpp +++ b/opm/core/props/satfunc/SatFuncStone2.hpp @@ -19,74 +19,180 @@ #ifndef SATFUNCSTONE2_HPP #define SATFUNCSTONE2_HPP -#include -#include -#include -#include -#include +#include namespace Opm { - class SatFuncStone2Uniform : public BlackoilPhases + template + class SatFuncStone2 : public SatFuncBase { public: - void init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples); void evalKr(const double* s, double* kr) const; void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalPc(const double* s, double* pc) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - double krwmax_; // Max water relperm - double kromax_; // Max oil relperm - double swcr_; // Critical water saturation. - double krwr_; // Water relperm at critical oil-in-water saturation. - double sowcr_; // Critical oil-in-water saturation. - double krorw_; // Oil relperm at critical water saturation. + + void evalKr(const double* s, double* kr, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} + void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const + {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const + {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} + void evalPc(const double* s, double* pc, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} + void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const + {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} + + private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - UniformTableLinear krw_; - UniformTableLinear krow_; - UniformTableLinear pcow_; - UniformTableLinear krg_; - UniformTableLinear krog_; - UniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) + }; + + typedef SatFuncStone2 > SatFuncStone2Uniform; + typedef SatFuncStone2 > SatFuncStone2Nonuniform; - - class SatFuncStone2Nonuniform : public BlackoilPhases + template + void SatFuncStone2::evalKr(const double* s, double* kr) const { - public: - void init(const EclipseGridParser& deck, - const int table_num, - const PhaseUsage phase_usg, - const int samples); - void evalKr(const double* s, double* kr) const; - void evalKrDeriv(const double* s, double* kr, double* dkrds) const; - void evalPc(const double* s, double* pc) const; - void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - double smin_[PhaseUsage::MaxNumPhases]; - double smax_[PhaseUsage::MaxNumPhases]; - double krwmax_; // Max water relperm - double kromax_; // Max oil relperm - double swcr_; // Critical water saturation. - double krwr_; // Water relperm at critical oil-in-water saturation. - double sowcr_; // Critical oil-in-water saturation. - double krorw_; // Oil relperm at critical water saturation. - private: - PhaseUsage phase_usage; // A copy of the outer class' phase_usage_. - NonuniformTableLinear krw_; - NonuniformTableLinear krow_; - NonuniformTableLinear pcow_; - NonuniformTableLinear krg_; - NonuniformTableLinear krog_; - NonuniformTableLinear pcog_; - double krocw_; // = krow_(s_wc) - }; + if (this->phase_usage.num_phases == 3) { + // Stone-II relative permeability model. + double sw = s[BlackoilPhases::Aqua]; + double sg = s[BlackoilPhases::Vapour]; + double krw = this->krw_(sw); + double krg = this->krg_(sg); + double krow = this->krow_(sw + sg); // = 1 - so + double krog = this->krog_(sg); // = 1 - so - sw + double krocw = this->krocw_; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); + if (kr[BlackoilPhases::Liquid] < 0.0) { + kr[BlackoilPhases::Liquid] = 0.0; + } + return; + } + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double krow = this->krow_(sw); + kr[wpos] = krw; + kr[opos] = krow; + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double krog = this->krog_(sg); + kr[gpos] = krg; + kr[opos] = krog; + } + } + + template + void SatFuncStone2::evalKrDeriv(const double* s, double* kr, double* dkrds) const + { + const int np = this->phase_usage.num_phases; + std::fill(dkrds, dkrds + np*np, 0.0); + + if (np == 3) { + // Stone-II relative permeability model. + double sw = s[BlackoilPhases::Aqua]; + double sg = s[BlackoilPhases::Vapour]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krow = this->krow_(sw + sg); + double dkrow = this->krow_.derivative(sw + sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + double krocw = this->krocw_; + kr[BlackoilPhases::Aqua] = krw; + kr[BlackoilPhases::Vapour] = krg; + kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg); + if (kr[BlackoilPhases::Liquid] < 0.0) { + kr[BlackoilPhases::Liquid] = 0.0; + } + dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww; + dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg; + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww); + dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg) + + krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg); + return; + } + // We have a two-phase situation. We know that oil is active. + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sw = s[wpos]; + double krw = this->krw_(sw); + double dkrww = this->krw_.derivative(sw); + double krow = this->krow_(sw); + double dkrow = this->krow_.derivative(sw); + kr[wpos] = krw; + kr[opos] = krow; + dkrds[wpos + wpos*np] = dkrww; + dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order. + } else { + assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]); + int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid]; + double sg = s[gpos]; + double krg = this->krg_(sg); + double dkrgg = this->krg_.derivative(sg); + double krog = this->krog_(sg); + double dkrog = this->krog_.derivative(sg); + kr[gpos] = krg; + kr[opos] = krog; + dkrds[gpos + gpos*np] = dkrgg; + dkrds[opos + gpos*np] = dkrog; + } + + } + + template + void SatFuncStone2::evalPc(const double* s, double* pc) const + { + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + pc[pos] = this->pcow_(s[pos]); + } + if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + pc[pos] = this->pcog_(s[pos]); + } + } + + template + void SatFuncStone2::evalPcDeriv(const double* s, double* pc, double* dpcds) const + { + // The problem of determining three-phase capillary pressures + // is very hard experimentally, usually one extends two-phase + // data (as for relative permeability). + // In our approach the derivative matrix is quite sparse, only + // the diagonal elements corresponding to non-oil phases are + // (potentially) nonzero. + const int np = this->phase_usage.num_phases; + std::fill(dpcds, dpcds + np*np, 0.0); + pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0; + if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua]; + pc[pos] = this->pcow_(s[pos]); + dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]); + } + if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) { + int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour]; + pc[pos] = this->pcog_(s[pos]); + dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]); + } + } } // namespace Opm #endif // SATFUNCSTONE2_HPP diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index f2892e12..c321a42a 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -101,34 +101,54 @@ 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. + void updateSatHyst(const int n, + const int* cells, + const double* s); private: PhaseUsage phase_usage_; std::vector satfuncset_; std::vector cell_to_func_; // = SATNUM - 1 - struct { // End point scaling parameters - std::vector swl_; - std::vector swcr_; - std::vector swu_; - std::vector sowcr_; - std::vector krw_; - std::vector krwr_; - std::vector kro_; - std::vector krorw_; - } eps_; bool do_eps_; // ENDSCALE is active bool do_3pt_; // SCALECRS: YES~true NO~false + bool do_hyst_; // Keywords ISWL etc detected + std::vector eps_transf_; + std::vector eps_transf_hyst_; + std::vector sat_hyst_; typedef SatFuncSet Funcs; const Funcs& funcForCell(const int cell) const; - void initEPS(const EclipseGridParser& deck, - const UnstructuredGrid& grid, - const std::string& keyword, - std::vector& scaleparam); - void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const; + const UnstructuredGrid& grid); + void initEPSHyst(const EclipseGridParser& deck, + const UnstructuredGrid& grid); + void initEPSKey(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam); + void initEPSParam(const int cell, + EPSTransforms::Transform& data, + const bool oil, + const double sl_tab, + const double scr_tab, + const double su_tab, + const double sxcr_tab, + const double s0_tab, + const double krsr_tab, + const double krmax_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); }; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index 0b5d6179..b2b9fd47 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -103,9 +103,9 @@ namespace Opm do_eps_ = false; do_3pt_ = false; if (deck.hasField("ENDSCALE")) { - if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) { - OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas."); - } + //if (!phase_usage_.phase_used[Aqua] || !phase_usage_.phase_used[Liquid] || phase_usage_.phase_used[Vapour]) { + // OPM_THROW(std::runtime_error, "Currently endpoint-scaling limited to oil-water systems without gas."); + //} if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) { OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted."); } @@ -118,14 +118,23 @@ namespace Opm } } do_eps_ = true; - initEPS(deck, grid, std::string("SWCR"), eps_.swcr_); - initEPS(deck, grid, std::string("SWL"), eps_.swl_); - initEPS(deck, grid, std::string("SWU"), eps_.swu_); - initEPS(deck, grid, std::string("SOWCR"), eps_.sowcr_); - initEPS(deck, grid, std::string("KRW"), eps_.krw_); - initEPS(deck, grid, std::string("KRWR"), eps_.krwr_); - initEPS(deck, grid, std::string("KRO"), eps_.kro_); - initEPS(deck, grid, std::string("KRORW"), eps_.krorw_); + + initEPS(deck, grid); + + // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR + do_hyst_ = deck.hasField("ISWL") || deck.hasField("ISWU") || deck.hasField("ISWCR") || deck.hasField("ISGL") || + deck.hasField("ISGU") || deck.hasField("ISGCR") || deck.hasField("ISOWCR") || deck.hasField("ISOGCR"); + if (do_hyst_) { + if (deck.hasField("KRW") || deck.hasField("KRG") || deck.hasField("KRO") || deck.hasField("KRWR") || + deck.hasField("KRGR") || deck.hasField("KRORW") || deck.hasField("KRORG") || + deck.hasField("IKRW") || deck.hasField("IKRG") || deck.hasField("IKRO") || deck.hasField("IKRWR") || + deck.hasField("IKRGR") || deck.hasField("IKRORW") || deck.hasField("IKRORG") ) { + OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Currently hysteresis and relperm value scaling can not be combined."); + } + initEPSHyst(deck, grid); + } + + //OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- ENDSCALE: Under construction ..."); } } @@ -165,8 +174,10 @@ namespace Opm if (dkrds) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - if (do_eps_) { - relpermEPS(s + np*i, cells[i], kr + np*i, dkrds + np*np*i); + if (do_hyst_) { + funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); + } else if (do_eps_) { + funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i, &(eps_transf_[cells[i]])); } else { funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); } @@ -174,8 +185,10 @@ namespace Opm } else { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - if (do_eps_) { - relpermEPS(s + np*i, cells[i], kr + np*i); + if (do_hyst_) { + funcForCell(cells[i]).evalKr(s + np*i, kr + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); + } else if (do_eps_) { + funcForCell(cells[i]).evalKr(s + np*i, kr + np*i, &(eps_transf_[cells[i]])); } else { funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); } @@ -209,12 +222,20 @@ namespace Opm if (dpcds) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i); + if (do_eps_) { + funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i, &(eps_transf_[cells[i]])); + } else { + funcForCell(cells[i]).evalPcDeriv(s + np*i, pc + np*i, dpcds + np*np*i); + } } } else { // #pragma omp parallel for - for (int i = 0; i < n; ++i) { - funcForCell(cells[i]).evalPc(s + np*i, pc + np*i); + for (int i = 0; i < n; ++i) { + if (do_eps_) { + funcForCell(cells[i]).evalPc(s + np*i, pc + np*i, &(eps_transf_[cells[i]])); + } else { + funcForCell(cells[i]).evalPc(s + np*i, pc + np*i); + } } } } @@ -234,16 +255,60 @@ namespace Opm double* smax) const { assert(cells != 0); - const int np = phase_usage_.num_phases; - for (int i = 0; i < n; ++i) { - for (int p = 0; p < np; ++p) { - smin[np*i + p] = funcForCell(cells[i]).smin_[p]; - smax[np*i + p] = funcForCell(cells[i]).smax_[p]; + + if (do_eps_) { + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + for (int i = 0; i < n; ++i) { + smin[np*i + opos] = 1.0; + smax[np*i + opos] = 1.0; + if (phase_usage_.phase_used[Aqua]) { + smin[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].wat.smin: funcForCell(cells[i]).smin_[wpos]; + smax[np*i + wpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].wat.smax: funcForCell(cells[i]).smax_[wpos]; + smin[np*i + opos] -= smax[np*i + wpos]; + smax[np*i + opos] -= smin[np*i + wpos]; + } + if (phase_usage_.phase_used[Vapour]) { + smin[np*i + gpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].gas.smin: funcForCell(cells[i]).smin_[gpos]; + smax[np*i + gpos] = eps_transf_[cells[i]].wat.doNotScale ? eps_transf_[cells[i]].gas.smax: funcForCell(cells[i]).smax_[gpos]; + smin[np*i + opos] -= smax[np*i + gpos]; + smax[np*i + opos] -= smin[np*i + gpos]; + } + if (phase_usage_.phase_used[Vapour] && phase_usage_.phase_used[Aqua]) { + smin[np*i + opos] = std::max(0.0,smin[np*i + opos]); + } + } + } else { + for (int i = 0; i < n; ++i) { + for (int p = 0; p < np; ++p) { + smin[np*i + p] = funcForCell(cells[i]).smin_[p]; + smax[np*i + p] = funcForCell(cells[i]).smax_[p]; + } } } } + + /// Update saturation state for the hysteresis tracking + /// \param[in] n Number of data points. + /// \param[in] s Array of nP saturation values. + template + void SaturationPropsFromDeck::updateSatHyst(const int n, + const int* cells, + const double* s) + { + assert(cells != 0); + + const int np = phase_usage_.num_phases; + if (do_hyst_) { +// #pragma omp parallel for + for (int i = 0; i < n; ++i) { + funcForCell(cells[i]).updateSatHyst(s + np*i, &(eps_transf_[cells[i]]), &(eps_transf_hyst_[cells[i]]), &(sat_hyst_[cells[i]])); + } + } + } // Map the cell number to the correct function set. template @@ -253,13 +318,147 @@ namespace Opm return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; } - // Initialize saturation scaling parameter + // Initialize saturation scaling parameters template void SaturationPropsFromDeck::initEPS(const EclipseGridParser& deck, - const UnstructuredGrid& grid, - const std::string& keyword, - std::vector& scaleparam) + const UnstructuredGrid& grid) { + std::vector swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr; + std::vector krw, krg, kro, krwr, krgr, krorw, krorg; + // Initialize saturation scaling parameter + initEPSKey(deck, grid, std::string("SWL"), swl); + initEPSKey(deck, grid, std::string("SWU"), swu); + initEPSKey(deck, grid, std::string("SWCR"), swcr); + initEPSKey(deck, grid, std::string("SGL"), sgl); + initEPSKey(deck, grid, std::string("SGU"), sgu); + initEPSKey(deck, grid, std::string("SGCR"), sgcr); + initEPSKey(deck, grid, std::string("SOWCR"), sowcr); + initEPSKey(deck, grid, std::string("SOGCR"), sogcr); + initEPSKey(deck, grid, std::string("KRW"), krw); + initEPSKey(deck, grid, std::string("KRG"), krg); + initEPSKey(deck, grid, std::string("KRO"), kro); + initEPSKey(deck, grid, std::string("KRWR"), krwr); + initEPSKey(deck, grid, std::string("KRGR"), krgr); + initEPSKey(deck, grid, std::string("KRORW"), krorw); + initEPSKey(deck, grid, std::string("KRORG"), krorg); + + eps_transf_.resize(grid.number_of_cells); + + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; + const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + + for (int cell = 0; cell < grid.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); + // ### 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); + } 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); + // ### 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); + } 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); + } + } + } + + + // Initialize hysteresis saturation scaling parameters + template + void SaturationPropsFromDeck::initEPSHyst(const EclipseGridParser& deck, + const UnstructuredGrid& grid) + { + std::vector iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr; + std::vector ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg; + // Initialize hysteresis saturation scaling parameters + initEPSKey(deck, grid, std::string("ISWL"), iswl); + initEPSKey(deck, grid, std::string("ISWU"), iswu); + initEPSKey(deck, grid, std::string("ISWCR"), iswcr); + initEPSKey(deck, grid, std::string("ISGL"), isgl); + initEPSKey(deck, grid, std::string("ISGU"), isgu); + initEPSKey(deck, grid, std::string("ISGCR"), isgcr); + initEPSKey(deck, grid, std::string("ISOWCR"), isowcr); + initEPSKey(deck, grid, std::string("ISOGCR"), isogcr); + initEPSKey(deck, grid, std::string("IKRW"), ikrw); + initEPSKey(deck, grid, std::string("IKRG"), ikrg); + initEPSKey(deck, grid, std::string("IKRO"), ikro); + initEPSKey(deck, grid, std::string("IKRWR"), ikrwr); + initEPSKey(deck, grid, std::string("IKRGR"), ikrgr); + initEPSKey(deck, grid, std::string("IKRORW"), ikrorw); + initEPSKey(deck, grid, std::string("IKRORG"), ikrorg); + + + eps_transf_hyst_.resize(grid.number_of_cells); + sat_hyst_.resize(grid.number_of_cells); + + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; + const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + + for (int cell = 0; cell < grid.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); + } + } + } + + // Initialize saturation scaling parameter + template + void SaturationPropsFromDeck::initEPSKey(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam) + { + const bool useAqua = phase_usage_.phase_used[Aqua]; + const bool useLiquid = phase_usage_.phase_used[Liquid]; + const bool useVapour = phase_usage_.phase_used[Vapour]; bool useKeyword = deck.hasField(keyword); bool hasENPTVD = deck.hasField("ENPTVD"); bool hasENKRVD = deck.hasField("ENKRVD"); @@ -269,70 +468,120 @@ namespace Opm // Active keyword assigned default values for each cell (in case of possible box-wise assignment) int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - if (keyword[0] == 'S' && (useKeyword || hasENPTVD)) { - if (keyword == std::string("SWL")) { - if (useKeyword || deck.getENPTVD().mask_[0]) { + int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) { + if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) { + if (useAqua && (useKeyword || deck.getENPTVD().mask_[0])) { itab = 1; scaleparam.resize(grid.number_of_cells); for (int i=0; i 0) { table = deck.getENPTVD().table_; } - } else if (keyword[0] == 'K' && (useKeyword || hasENKRVD)) { - if (keyword == std::string("KRW")) { - if (useKeyword || deck.getENKRVD().mask_[0]) { + } else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) { + if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) { + if (useAqua && (useKeyword || deck.getENKRVD().mask_[0])) { itab = 1; scaleparam.resize(grid.number_of_cells); for (int i=0; i - void SaturationPropsFromDeck::relpermEPS(const double *s, const int cell, double *kr, double *dkrds) const + void SaturationPropsFromDeck::initEPSParam(const int cell, + EPSTransforms::Transform& data, + const bool oil, // flag indicating krow/krog calculations + const double sl_tab, // minimum saturation (for krow/krog calculations this is normally zero) + const double scr_tab, // critical saturation + const double su_tab, // maximum saturation (for krow/krog calculations this is minimum water/gas saturation) + const double sxcr_tab, // second critical saturation (not used for 2pt scaling) + const double s0_tab, // threephase complementary minimum saturation (-1.0 indicates 2-phase) + const double krsr_tab, // relperm at displacing critical saturation + const double krmax_tab, // relperm at maximum saturation + 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 int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; - const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; - double ss[PhaseUsage::MaxNumPhases]; + if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) { + data.doNotScale = true; + } else { + data.doNotScale = false; + data.do_3pt = do_3pt_; + double s_r; + if (s0_tab < 0.0) { // 2phase + s_r = 1.0-sxcr_tab; + if (do_3pt_) data.sr = sxcr.empty() ? s_r : 1.0-sxcr[cell]; + } else { // 3phase + s_r = 1.0-sxcr_tab-s0_tab; + if (do_3pt_)data.sr = 1.0 - (sxcr.empty() ? sxcr_tab : sxcr[cell]) + - (s0.empty() ? s0_tab : s0[cell]); + } + data.scr = scr.empty() ? scr_tab : scr[cell]; + double s_max = su_tab; + if (oil) { + data.smin = sl_tab; + if (s0_tab < 0.0) { // 2phase + s_max = 1.0 - su_tab; + data.smax = 1.0 - (su.empty() ? su_tab : su[cell]); + } else { // 3phase + s_max = 1.0 - su_tab - s0_tab; + data.smax = 1.0 - (su.empty() ? su_tab : su[cell]) + - (s0.empty() ? s0_tab : s0[cell]); + } + } else { + data.smin = sl.empty() ? sl_tab : sl[cell]; + data.smax = su.empty() ? su_tab : su[cell]; + } + if (do_3pt_) { + data.slope1 = (s_r-scr_tab)/(data.sr-data.scr); + data.slope2 = (s_max-s_r)/(data.smax-data.sr); + } else { + data.slope2 = data.slope1 = (s_max-scr_tab)/(data.smax-data.scr); + // Inv transform of tabulated critical displacing saturation to prepare for possible value scaling (krwr etc) + data.sr = data.scr + (s_r-scr_tab)*(data.smax-data.scr)/(s_max-scr_tab); + } + } + + data.doKrMax = !krmax.empty(); + data.doKrCrit = !krsr.empty(); + data.doSatInterp = false; + data.krsr = krsr.empty() ? krsr_tab : krsr[cell]; + data.krmax = krmax.empty() ? krmax_tab : krmax[cell]; + data.krSlopeCrit = data.krsr/krsr_tab; + data.krSlopeMax = data.krmax/krmax_tab; + if (data.doKrCrit) { + if (data.sr > data.smax-1.0e-6) { + //Ignore krsr and do two-point (one might consider combining krsr and krmax linearly between scr and smax ... ) + data.doKrCrit = false; + } else if (std::fabs(krmax_tab- krsr_tab) > 1.0e-6) { // interpolate kr + data.krSlopeMax = (data.krmax-data.krsr)/(krmax_tab-krsr_tab); + } else { // interpolate sat + data.doSatInterp = true; + data.krSlopeMax = (data.krmax-data.krsr)/(data.smax-data.sr); + } + } - if (do_3pt_) { // Three-point scaling - // Transforms for water saturation - if (eps_.swcr_.empty() && eps_.swu_.empty()) { - ss[wpos] = s[wpos]; - } else { - double s_r = 1.0-funcForCell(cell).sowcr_; - double sr = eps_.sowcr_.empty() ? s_r : 1.0-eps_.sowcr_[cell]; - if (s[wpos] <= sr) { - double sw_cr = funcForCell(cell).swcr_; - double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell]; - ss[wpos] = (s[wpos] <= swcr) ? sw_cr : sw_cr+(s[wpos]-swcr)*(s_r-sw_cr)/(sr-swcr); - } else { - double sw_max = funcForCell(cell).smax_[wpos]; - double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell]; - ss[wpos] = (s[wpos] >= swmax) ? sw_max : s_r+(s[wpos]-sr)*(sw_max-s_r)/(swmax-sr); - } - } - // Transforms for oil saturation - if (eps_.sowcr_.empty() && eps_.swl_.empty()) { - ss[opos] = s[opos]; - } else { - double s_r = 1.0-funcForCell(cell).swcr_; - double sr = eps_.swcr_.empty() ? s_r : 1.0-eps_.swcr_[cell]; - if (s[opos] <= sr) { - double sow_cr = funcForCell(cell).sowcr_; - double sowcr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell]; - ss[opos] = (s[opos] <= sowcr) ? sow_cr : sow_cr+(s[opos]-sowcr)*(s_r-sow_cr)/(sr-sowcr); - } else { - double sow_max = funcForCell(cell).smax_[opos]; - double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]); - ss[opos] = (s[opos] >= sowmax) ? sow_max : s_r+(s[opos]-sr)*(sow_max-s_r)/(sowmax-sr); - } - } - } else { // Two-point scaling - // Transforms for water saturation - if (eps_.swcr_.empty() && eps_.swu_.empty()) { - ss[wpos] = s[wpos]; - } else { - double sw_cr = funcForCell(cell).swcr_; - double swcr = eps_.swcr_.empty() ? sw_cr : eps_.swcr_[cell]; - if (s[wpos] <= swcr) { - ss[wpos] = sw_cr; - } else { - double sw_max = funcForCell(cell).smax_[wpos]; - double swmax = eps_.swu_.empty() ? sw_max : eps_.swu_[cell]; - ss[wpos] = (s[wpos] >= swmax) ? sw_max : sw_cr + (s[wpos]-swcr)*(sw_max-sw_cr)/(swmax-swcr); - } - } - // Transforms for oil saturation - if (eps_.sowcr_.empty() && eps_.swl_.empty()) { - ss[opos] = s[opos]; - } else { - double sow_cr = funcForCell(cell).sowcr_; - double socr = eps_.sowcr_.empty() ? sow_cr : eps_.sowcr_[cell]; - if (s[opos] <= socr) { - ss[opos] = sow_cr; - } else { - double sow_max = funcForCell(cell).smax_[opos]; - double sowmax = eps_.swl_.empty() ? sow_max : (1.0-eps_.swl_[cell]); - ss[opos] = (s[opos] >= sowmax) ? sow_max : sow_cr + (s[opos]-socr) *(sow_max-sow_cr)/(sowmax-socr); - } - } - } - - // Evaluation of relperms - if (dkrds) { - OPM_THROW(std::runtime_error, "Relperm derivatives not yet available in combination with end point scaling ..."); - funcForCell(cell).evalKrDeriv(ss, kr, dkrds); - } else { - // Assume: sw_cr -> krw=0 sw_max -> krw= - // sow_cr -> kro=0 sow_max -> kro= - funcForCell(cell).evalKr(ss, kr); - } - - // Scaling of relperms values - // - Water - if (eps_.krw_.empty() && eps_.krwr_.empty()) { // No value scaling - } else if (eps_.krwr_.empty()) { // Two-point - kr[wpos] *= (eps_.krw_[cell]/funcForCell(cell).krwmax_); - } else { - double swcr = eps_.swcr_.empty() ? funcForCell(cell).swcr_ : eps_.swcr_[cell]; - double swmax = eps_.swu_.empty() ? funcForCell(cell).smax_[wpos] : eps_.swu_[cell]; - double sr; - if (do_3pt_) { - sr = eps_.sowcr_.empty() ? 1.0-funcForCell(cell).sowcr_ : 1.0-eps_.sowcr_[cell]; - } else { - double sw_cr = funcForCell(cell).swcr_; - double sw_max = funcForCell(cell).smax_[wpos]; - double s_r = 1.0-funcForCell(cell).sowcr_; - sr = swcr + (s_r-sw_cr)*(swmax-swcr)/(sw_max-sw_cr); - } - if (s[wpos] <= swcr) { - kr[wpos] = 0.0; - } else if (sr > swmax-1.0e-6) { - if (do_3pt_) { //Ignore krw and do two-point? - kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_; - } else if (!eps_.kro_.empty()){ //Ignore krwr and do two-point - kr[wpos] *= eps_.krw_[cell]/funcForCell(cell).krwmax_; - } - } else if (s[wpos] <= sr) { - kr[wpos] *= eps_.krwr_[cell]/funcForCell(cell).krwr_; - } else if (s[wpos] <= swmax) { - double krw_max = funcForCell(cell).krwmax_; - double krw = eps_.krw_.empty() ? krw_max : eps_.krw_[cell]; - double krw_r = funcForCell(cell).krwr_; - double krwr = eps_.krwr_.empty() ? krw_r : eps_.krwr_[cell]; - if (std::fabs(krw_max- krw_r) > 1.0e-6) { - kr[wpos] = krwr + (kr[wpos]-krw_r)*(krw-krwr)/(krw_max-krw_r); - } else { - kr[wpos] = krwr + (krw-krwr)*(s[wpos]-sr)/(swmax-sr); - } - } else { - kr[wpos] = eps_.krw_.empty() ? funcForCell(cell).krwmax_ : eps_.krw_[cell]; - } - } - - // - Oil - if (eps_.kro_.empty() && eps_.krorw_.empty()) { // No value scaling - } else if (eps_.krorw_.empty()) { // Two-point scaling - kr[opos] *= (eps_.kro_[cell]/funcForCell(cell).kromax_); - } else { - double sowcr = eps_.sowcr_.empty() ? funcForCell(cell).sowcr_ : eps_.sowcr_[cell]; - double sowmax = eps_.swl_.empty() ? funcForCell(cell).smax_[opos] : 1.0-eps_.swl_[cell]; - double sr; - if (do_3pt_) { - sr = eps_.swcr_.empty() ? 1.0-funcForCell(cell).swcr_ : 1.0-eps_.swcr_[cell]; - } else { - double sow_cr = funcForCell(cell).sowcr_; - double sow_max = funcForCell(cell).smax_[opos]; - double s_r = 1.0-funcForCell(cell).swcr_; - sr = sowcr + (s_r-sow_cr)*(sowmax-sowcr)/(sow_max-sow_cr); - } - if (s[opos] <= sowcr) { - kr[opos] = 0.0; - } else if (sr > sowmax-1.0e-6) { - if (do_3pt_) { //Ignore kro and do two-point? - kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_; - } else if (!eps_.kro_.empty()){ //Ignore krowr and do two-point - kr[opos] *= eps_.kro_[cell]/funcForCell(cell).kromax_; - } - } else if (s[opos] <= sr) { - kr[opos] *= eps_.krorw_[cell]/funcForCell(cell).krorw_; - } else if (s[opos] <= sowmax) { - double kro_max = funcForCell(cell).kromax_; - double kro = eps_.kro_.empty() ? kro_max : eps_.kro_[cell]; - double kro_rw = funcForCell(cell).krorw_; - double krorw = eps_.krorw_[cell]; - if (std::fabs(kro_max- kro_rw) > 1.0e-6) { - kr[opos] = krorw + (kr[opos]- kro_rw)*(kro-krorw)/(kro_max- kro_rw); - } else { - kr[opos] = krorw + (kro-krorw)*(s[opos]-sr)/(sowmax-sr); - } - } else { - kr[opos] = eps_.kro_.empty() ? funcForCell(cell).kromax_ : eps_.kro_[cell]; - } - } } } // namespace Opm