diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f9d502c0..963a9a2a 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -206,6 +206,7 @@ list (APPEND TEST_DATA_FILES tests/satfuncEPS_A.DATA tests/satfuncEPS_B.DATA tests/satfuncEPS_C.DATA + tests/satfuncEPS_D.DATA tests/testBlackoilState1.DATA tests/testBlackoilState2.DATA tests/wells_manager_data.data diff --git a/opm/core/props/satfunc/SatFuncGwseg.hpp b/opm/core/props/satfunc/SatFuncGwseg.hpp index b35fefdf..4a6d0fed 100644 --- a/opm/core/props/satfunc/SatFuncGwseg.hpp +++ b/opm/core/props/satfunc/SatFuncGwseg.hpp @@ -286,7 +286,7 @@ namespace Opm 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); + 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; @@ -469,21 +469,21 @@ namespace Opm 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)); + krow = epst->watoil.scaleKr(ssow, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_); + dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(std::max(1.0-_ssow,_swco))); } 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)); + krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(std::max(1.0-_ssow,_swco)), this->krorw_); + dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(std::max(1.0-_ssow,_swco))); } // 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)); + double krog = epst->gasoil.scaleKr(ssog, this->krog_(std::max(1.0-_ssog-_swco,eps)), this->krorg_); + double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(std::max(1.0-_ssog-_swco,eps))); // xw and xg are the fractions occupied by water and gas zones. const double xw = (sw - swco) / d; diff --git a/tests/satfuncEPS_D.DATA b/tests/satfuncEPS_D.DATA new file mode 100644 index 00000000..1077121e --- /dev/null +++ b/tests/satfuncEPS_D.DATA @@ -0,0 +1,181 @@ +NOECHO + +RUNSPEC ====== + +WATER +OIL +GAS +DISGAS +VAPOIL + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 10 +/ + +WELLDIMS + 30 10 2 30 / + +ENDSCALE +--DIR REV NTENDP NSENDP +'NODIR' 'REVERS' 1 20 / +/ + +SATOPTS +HYSTER / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +10*5.0 +/ + + +PORO +10*0.2 +/ + + +PERMZ + 10*1.0 +/ + +PERMY +10*100.0 +/ + +PERMX +10*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + +-- Hysteresis input +EHYSTR + 0.1 0 0.1 1* KR / + +PVTO +-- Rs Pbub Bo Vo + 0 1. 1.0000 1.20 / + 20 40. 1.0120 1.17 / + 40 80. 1.0255 1.14 / + 60 120. 1.0380 1.11 / + 80 160. 1.0510 1.08 / + 100 200. 1.0630 1.06 / + 120 240. 1.0750 1.03 / + 140 280. 1.0870 1.00 / + 160 320. 1.0985 .98 / + 180 360. 1.1100 .95 / + 200 400. 1.1200 .94 + 500. 1.1189 .94 / +/ + +PVTG +-- Pg Rv Bg Vg + 100 0.0001 0.010 0.1 + 0.0 0.0104 0.1 / + 200 0.0004 0.005 0.2 + 0.0 0.0054 0.2 / +/ + +SCALECRS +-- YES / + NO/ + +SWOF +0.1 0.0 1.0 0.9 +0.2 0.0 0.8 0.8 +0.3 0.1 0.6 0.7 +0.4 0.2 0.4 0.6 +0.7 0.5 0.1 0.3 +0.8 0.6 0.0 0.2 +0.9 0.7 0.0 0.1 +/ + +SGOF +0.0 0.0 1.0 0.2 +0.1 0.0 0.7 0.4 +0.2 0.1 0.6 0.6 +0.8 0.7 0.0 2.0 +0.9 1.0 0.0 2.1 +/ + +PVTW +--RefPres Bw Comp Vw Cv + 1. 1.0 4.0E-5 0.96 0.0 / + + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +DENSITY +700 1000 1 +/ + +REGIONS ====== + +ENDNUM +10*1 +/ + +SOLUTION ====== + +EQUIL +45 150 50 0.25 45 0.35 1 1 0 +/ + +RSVD + 0 0.0 + 100 100. / + +RVVD + 0. 0. + 100. 0.0001 / + +RPTSOL +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +RPTSCHED +'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/test_satfunc.cpp b/tests/test_satfunc.cpp index 15ebac46..406ca81f 100644 --- a/tests/test_satfunc.cpp +++ b/tests/test_satfunc.cpp @@ -219,7 +219,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPSBase) BOOST_AUTO_TEST_CASE (GwsegEPS_A) { // This is the eps (but no hysteris) version of the Gwseg model. - // However, only default scaling parameters, i.e no scaling. + // Scaling parameters from keyword SWL etc. //std::cout << "==================================== GwsegEPS_A ====================================" << std::endl; @@ -366,7 +366,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_A) BOOST_AUTO_TEST_CASE (GwsegEPS_B) { // This is the eps (but no hysteris) version of the Gwseg model. - // However, only default scaling parameters, i.e no scaling. + // Scaling parameters from ENPTVD table. //std::cout << "==================================== GwsegEPS_B ====================================" << std::endl; /* @@ -465,7 +465,7 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_B) BOOST_AUTO_TEST_CASE (GwsegEPS_C) { // This is the eps (but no hysteris) version of the Gwseg model. - // However, only default scaling parameters, i.e no scaling. + // Scaling parameters given the "norne-way", i.e EQUALS, COPY, ADD, MULTIPLY. //std::cout << "==================================== GwsegEPS_C ====================================" << std::endl; @@ -560,4 +560,89 @@ BOOST_AUTO_TEST_CASE (GwsegEPS_C) } +BOOST_AUTO_TEST_CASE (GwsegEPS_D) +{ + // This is the eps and hysteris version of the Gwseg model. + // However, only default scaling parameters, i.e no scaling. + + //std::cout << "==================================== GwsegEPS_D ====================================" << std::endl; + + Opm::parameter::ParameterGroup param; + + Opm::GridManager gm(1, 1, 10, 1.0, 1.0, 5.0); + const UnstructuredGrid& grid = *(gm.c_grid()); + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile("satfuncEPS_D.DATA"); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(deck)); + Opm::BlackoilPropertiesFromDeck props(deck, eclipseState, grid, param, false); + + const int np = props.numPhases(); + const int wpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Aqua]; + const int opos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Liquid]; + const int gpos = props.phaseUsage().phase_pos[Opm::BlackoilPhases::Vapour]; + + BOOST_REQUIRE(np == 3); + BOOST_REQUIRE(wpos == 0); + BOOST_REQUIRE(opos == 1); + BOOST_REQUIRE(gpos == 2); + + const int n=11; + double s[n*np]; + int cells[n]; + double kr[n*np]; + double dkrds[n*np*np]; + + for (int i=0; i