diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index 9c0f437d..fef38c44 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -97,7 +97,10 @@ namespace EclipseKeywords string("MULTPV"), string("PRESSURE"), string("SGAS"), string("SWAT"), string("SOIL"), string("RS"), string("DXV"), string("DYV"), string("DZV"), - string("DEPTHZ"), string("TOPS"), string("MAPAXES") + string("DEPTHZ"), string("TOPS"), string("MAPAXES"), + string("SWCR"), string("SWL"), string("SWU"), + string("SOWCR"), string("KRW"), string("KRWR"), + string("KRO"), string("KRORW") }; const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); @@ -114,7 +117,7 @@ namespace EclipseKeywords string("PLYVISC"), string("PLYROCK"), string("PLYADS"), string("PLYMAX"), string("TLMIXPAR"), string("WPOLYMER"), string("GRUPTREE"), string("GCONINJE"), string("GCONPROD"), - string("WGRUPCON"), + string("WGRUPCON"), string("ENDSCALE"), string("SCALECRS"), // The following fields only have a dummy implementation // that allows us to ignore them. string("SWFN"), @@ -554,7 +557,9 @@ void EclipseGridParser::convertToSI() key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" || key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" || key == "SGAS" || key == "SWAT" || key == "SOIL" || - key == "RS") { + key == "RS" || key == "SWCR" || key == "SWL" || + key == "SWU" || key == "SOWCR" || key == "KRW" || + key == "KRWR" || key == "KRORW" || key == "KRO") { unit = 1.0; do_convert = false; // Dimensionless keywords... } else if (key == "PRESSURE") { diff --git a/opm/core/eclipse/EclipseGridParser.hpp b/opm/core/eclipse/EclipseGridParser.hpp index e38f8b28..ab871c68 100644 --- a/opm/core/eclipse/EclipseGridParser.hpp +++ b/opm/core/eclipse/EclipseGridParser.hpp @@ -192,6 +192,8 @@ public: SPECIAL_FIELD(GCONINJE) SPECIAL_FIELD(GCONPROD) SPECIAL_FIELD(WGRUPCON) + SPECIAL_FIELD(ENDSCALE) + SPECIAL_FIELD(SCALECRS) // The following fields only have a dummy implementation // that allows us to ignore them. diff --git a/opm/core/eclipse/SpecialEclipseFields.hpp b/opm/core/eclipse/SpecialEclipseFields.hpp index 507073d3..70bff0a0 100644 --- a/opm/core/eclipse/SpecialEclipseFields.hpp +++ b/opm/core/eclipse/SpecialEclipseFields.hpp @@ -2206,6 +2206,109 @@ struct WPOLYMER : public SpecialBase } }; +struct ENDSCALE : public SpecialBase { + std::string dir_switch_; + std::string revers_switch_; + int n_tab_; + int n_node_; + + virtual std::string name() const { + return std::string("ENDSCALE"); + } + + virtual void read(std::istream & is) { + dir_switch_ = std::string("NODIR"); + revers_switch_ = std::string("REVERS"); + n_tab_ = 1; + n_node_ = 20; + while (!is.eof()) { + std::string tmp = readString(is); + while (tmp.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + tmp = readString(is); + } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } + dir_switch_ = tmp; + + tmp = readString(is); + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } + revers_switch_ = tmp; + + is >> ignoreWhitespace; + if(is.peek() == int('/')) { + is >> ignoreLine; + return; + } + is >> n_tab_; + + is >> ignoreWhitespace; + if(is.peek() == int('/')) { + is >> ignoreLine; + return; + } + is >> n_node_; + + is >> ignoreLine; + break; + } + } + + virtual void write(std::ostream & os) { + os << name() << '\n'; + os << dir_switch_ << " " << revers_switch_ << " " + << n_tab_ << " " << n_node_ << '\n'; + } + + virtual void convertToSI(const EclipseUnits&) { + } +}; + + +struct SCALECRS : public SpecialBase { + std::string scalecrs_; + + virtual std::string name() const { + return std::string("SCALECRS"); + } + + virtual void read(std::istream & is) { + scalecrs_ = std::string("NO"); + while (!is.eof()) { + std::string tmp = readString(is); + while (tmp.find("--") == 0) { + // This line is a comment + is >> ignoreLine; + tmp = readString(is); + } + if (tmp[0] == '/') { + is >> ignoreLine; + break; + } else if (tmp[0] == 'Y') { + scalecrs_ = std::string("YES"); + } + + is >> ignoreLine; + break; + } + } + + virtual void write(std::ostream & os) { + os << name() << '\n'; + os << scalecrs_ << '\n'; + } + + virtual void convertToSI(const EclipseUnits&) { + } +}; + + // The following fields only have a dummy implementation // that allows us to ignore them. struct SWFN : public MultRec {}; diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index 271f9395..d49dc5b9 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -50,6 +50,9 @@ namespace Opm // Unfortunate lack of pointer smartness here... const int sat_samples = param.getDefault("sat_tab_size", 200); std::string threephase_model = param.getDefault("threephase_model", "simple"); + if (deck.hasField("ENDSCALE") && threephase_model != "simple") { + THROW("Sorry, end point scaling currently available for the 'simple' model only."); + } if (sat_samples > 1) { if (threephase_model == "stone2") { SaturationPropsFromDeck* ptr diff --git a/opm/core/fluid/SatFuncGwseg.hpp b/opm/core/fluid/SatFuncGwseg.hpp index 1932e7f0..ae5b47a6 100644 --- a/opm/core/fluid/SatFuncGwseg.hpp +++ b/opm/core/fluid/SatFuncGwseg.hpp @@ -40,6 +40,12 @@ namespace Opm 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_. UniformTableLinear krw_; @@ -65,6 +71,12 @@ namespace Opm 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_; diff --git a/opm/core/fluid/SatFuncSimple.cpp b/opm/core/fluid/SatFuncSimple.cpp index 7e226a4a..37e74fff 100644 --- a/opm/core/fluid/SatFuncSimple.cpp +++ b/opm/core/fluid/SatFuncSimple.cpp @@ -54,6 +54,27 @@ namespace Opm 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 (unsigned int i=1; i 0.0) { + swcr_ = sw[i-1]; + krorw_ = krow[i-1]; + break; + } + } + for (unsigned int i=sw.size()-2; i>=0; --i) { + if (krow[i]> 0.0) { + sowcr_ = 1.0 - sw[i+1]; + krwr_ = krw[i+1]; + break; + } + } } if (phase_usage.phase_used[Vapour]) { const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; @@ -102,7 +123,8 @@ namespace Opm int opos = phase_usage.phase_pos[Liquid]; double sw = s[wpos]; double krw = krw_(sw); - double krow = krow_(sw); + double so = s[opos]; + double krow = krow_(1.0-so); kr[wpos] = krw; kr[opos] = krow; } else { @@ -160,8 +182,9 @@ namespace Opm double sw = s[wpos]; double krw = krw_(sw); double dkrww = krw_.derivative(sw); - double krow = krow_(sw); - double dkrow = krow_.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; @@ -253,6 +276,28 @@ namespace Opm 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 (unsigned int i=1; i 0.0) { + swcr_ = sw[i-1]; + krorw_ = krow[i-1]; + break; + } + } + for (unsigned int i=sw.size()-2; i>=0; --i) { + if (krow[i]> 0.0) { + sowcr_ = 1.0 - sw[i+1]; + krwr_ = krw[i+1]; + break; + } + } + } if (phase_usage.phase_used[Vapour]) { const SGOF::table_t& sgof_table = deck.getSGOF().sgof_; @@ -301,7 +346,8 @@ namespace Opm int opos = phase_usage.phase_pos[Liquid]; double sw = s[wpos]; double krw = krw_(sw); - double krow = krow_(sw); + double so = s[opos]; + double krow = krow_(1.0-so); kr[wpos] = krw; kr[opos] = krow; } else { @@ -359,8 +405,9 @@ namespace Opm double sw = s[wpos]; double krw = krw_(sw); double dkrww = krw_.derivative(sw); - double krow = krow_(sw); - double dkrow = krow_.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; diff --git a/opm/core/fluid/SatFuncSimple.hpp b/opm/core/fluid/SatFuncSimple.hpp index 338a359e..1b6ba03b 100644 --- a/opm/core/fluid/SatFuncSimple.hpp +++ b/opm/core/fluid/SatFuncSimple.hpp @@ -40,6 +40,12 @@ namespace Opm 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_. UniformTableLinear krw_; @@ -65,6 +71,12 @@ namespace Opm 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_; diff --git a/opm/core/fluid/SatFuncStone2.hpp b/opm/core/fluid/SatFuncStone2.hpp index 532bf7ce..fc906126 100644 --- a/opm/core/fluid/SatFuncStone2.hpp +++ b/opm/core/fluid/SatFuncStone2.hpp @@ -40,6 +40,12 @@ namespace Opm 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_. UniformTableLinear krw_; @@ -65,6 +71,12 @@ namespace Opm 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_; diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index ea2acb34..1492c5a5 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -107,9 +107,28 @@ namespace Opm 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 + 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; }; diff --git a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp index f7a80b45..6b87bc29 100644 --- a/opm/core/fluid/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck_impl.hpp @@ -96,6 +96,61 @@ namespace Opm for (int table = 0; table < num_tables; ++table) { satfuncset_[table].init(deck, table, phase_usage_, samples); } + + // Saturation table scaling + 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]) { + THROW("Currently endpoint-scaling limited to oil-water systems without gas."); + } + if (deck.getENDSCALE().dir_switch_ != std::string("NODIR")) { + THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'NODIR' accepted."); + } + if (deck.getENDSCALE().revers_switch_ != std::string("REVERS")) { + THROW("SaturationPropsFromDeck::init() -- ENDSCALE: Currently only 'REVERS' accepted."); + } + if (deck.hasField("SCALECRS")) { + if (deck.getSCALECRS().scalecrs_ == std::string("YES")) { + do_3pt_ = true; + } + } + 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_); + +/* + double ss[PhaseUsage::MaxNumPhases], kr[PhaseUsage::MaxNumPhases]; + int oldP = std::cout.precision(); + std::cout.precision(4); + for (unsigned int i=0; i<=100; ++i) { + ss[phase_usage_.phase_pos[Aqua]] = i*0.01; + ss[phase_usage_.phase_pos[Liquid]] = 1.0 - ss[phase_usage_.phase_pos[Aqua]]; + endScaling(ss, 15, kr); + std::cout << std::showpoint + << std::setw(10) << ss[phase_usage_.phase_pos[Aqua]] + << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]; + endScaling(ss, 45, kr); + std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]; + endScaling(ss, 75, kr); + std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]]; + endScaling(ss, 105, kr); + std::cout << std::setw(12) << kr[phase_usage_.phase_pos[Aqua]] + << std::setw(10) << kr[phase_usage_.phase_pos[Liquid]] + << std::noshowpoint << std::endl; + } + std::cout.precision(oldP); +*/ + } } @@ -134,12 +189,20 @@ namespace Opm if (dkrds) { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); + if (do_eps_) { + relpermEPS(s + np*i, cells[i], kr + np*i, dkrds + np*np*i); + } else { + funcForCell(cells[i]).evalKrDeriv(s + np*i, kr + np*i, dkrds + np*np*i); + } } } else { // #pragma omp parallel for for (int i = 0; i < n; ++i) { - funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); + if (do_eps_) { + relpermEPS(s + np*i, cells[i], kr + np*i); + } else { + funcForCell(cells[i]).evalKr(s + np*i, kr + np*i); + } } } } @@ -214,7 +277,223 @@ namespace Opm return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]]; } + // Initialize saturation scaling parameter + template + void SaturationPropsFromDeck::initEPS(const EclipseGridParser& deck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam) + { + + if (deck.hasField(keyword)) { + // Active keyword assigned default values for each cell (in case of possible box-wise assignment) + scaleparam.resize(grid.number_of_cells); + int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + if (keyword == std::string("SWCR")) { + for (int i=0; i& val = deck.getFloatingPointValue(keyword); + for (int c = 0; c < int(scaleparam.size()); ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + scaleparam[c] = val[deck_pos]; + } + } + } + + // Saturation scaling + template + void SaturationPropsFromDeck::relpermEPS(const double *s, const int cell, double *kr, double *dkrds) const + { + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int opos = phase_usage_.phase_pos[BlackoilPhases::Liquid]; + double ss[PhaseUsage::MaxNumPhases]; + + 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) { + THROW("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