diff --git a/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.hpp b/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.hpp index e81243978..66082785a 100644 --- a/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.hpp +++ b/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.hpp @@ -19,6 +19,7 @@ #ifndef ECLIPSE_SATFUNCPROPERTY_INITIALIZERS_HPP #define ECLIPSE_SATFUNCPROPERTY_INITIALIZERS_HPP +#include #include #include @@ -29,9 +30,35 @@ namespace Opm { namespace Opm { namespace satfunc { + struct RawTableEndPoints + { + struct { + std::vector gas; + std::vector water; + } connate; + + struct { + std::vector oil_in_gas; + std::vector oil_in_water; + std::vector gas; + std::vector water; + } critical; + + struct { + std::vector gas; + std::vector water; + } maximum; + }; + + std::shared_ptr + getRawTableEndpoints(const Opm::TableManager& tm, + const Opm::Phases& phases, + const double tolcrit); + std::vector init(const std::string& kewyord, const TableManager& tables, const Phases& phases, + const RawTableEndPoints& ep, const std::vector& cell_depth, const std::vector& num, const std::vector& endnum); diff --git a/src/opm/output/eclipse/Tables.cpp b/src/opm/output/eclipse/Tables.cpp index dca0d4633..b83c6c497 100644 --- a/src/opm/output/eclipse/Tables.cpp +++ b/src/opm/output/eclipse/Tables.cpp @@ -232,7 +232,7 @@ namespace { namespace SatFunc { std::transform(begin, end, dest, [tolcrit](const double kr) -> double { - return (kr < tolcrit) ? 0.0 : kr; + return (kr > tolcrit) ? kr : 0.0; }); } diff --git a/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.cpp b/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.cpp index 201a27c5a..9e3665b9b 100644 --- a/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.cpp +++ b/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.cpp @@ -409,6 +409,7 @@ FieldProps::FieldProps(const Deck& deck, const Phases& phases, const EclipseGrid ny(grid.getNY()), nz(grid.getNZ()), m_phases(phases), + m_satfuncctrl(deck), m_actnum(grid.getACTNUM()), cell_volume(extract_cell_volume(grid)), cell_depth(extract_cell_depth(grid)), @@ -1056,14 +1057,16 @@ void FieldProps::scanEDITSection(const EDITSection& edit_section) { void FieldProps::init_satfunc(const std::string& keyword, FieldData& satfunc) { + if (this->m_rtep == nullptr) + this->m_rtep = satfunc::getRawTableEndpoints(this->tables, this->m_phases, + this->m_satfuncctrl.minimumRelpermMobilityThreshold()); + const auto& endnum = this->get("ENDNUM"); - if (keyword[0] == 'I') { - const auto& imbnum = this->get("IMBNUM"); - satfunc.default_update(satfunc::init(keyword, this->tables, this->m_phases, this->cell_depth, imbnum, endnum)); - } else { - const auto& satnum = this->get("SATNUM"); - satfunc.default_update(satfunc::init(keyword, this->tables, this->m_phases, this->cell_depth, satnum, endnum)); - } + const auto& satreg = (keyword[0] == 'I') + ? this->get("IMBNUM") + : this->get("SATNUM"); + + satfunc.default_update(satfunc::init(keyword, this->tables, this->m_phases, *this->m_rtep, this->cell_depth, satreg, endnum)); } diff --git a/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.hpp b/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.hpp index fc9aac776..31d474c9f 100644 --- a/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.hpp +++ b/src/opm/parser/eclipse/EclipseState/Grid/FieldProps.hpp @@ -19,6 +19,7 @@ #ifndef FIELDPROPS_HPP #define FIELDPROPS_HPP +#include #include #include #include @@ -27,6 +28,7 @@ #include #include #include +#include #include namespace Opm { @@ -357,12 +359,14 @@ private: const UnitSystem unit_system; std::size_t nx,ny,nz; Phases m_phases; + SatFuncControls m_satfuncctrl; std::vector m_actnum; std::vector cell_volume; std::vector cell_depth; const std::string m_default_region; const EclipseGrid * grid_ptr; // A bit undecided whether to properly use the grid or not ... const TableManager& tables; + std::shared_ptr m_rtep; std::vector multregp; std::unordered_map> int_data; std::unordered_map> double_data; diff --git a/src/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.cpp b/src/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.cpp index 9bdc639b6..d117c986e 100644 --- a/src/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.cpp +++ b/src/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.cpp @@ -27,6 +27,7 @@ #include #include #include +#include #include #include @@ -34,34 +35,32 @@ #include #include +#include #include +#include #include +#include #include #include +#include +#include #include +// Note on deriving critical saturations: All table scanners are implemented +// in terms of std::lower_bound(begin, end, tolcrit, predicate) which returns +// the first position in [begin, end) for which +// +// predicate(*iter, tolcrit) +// +// is false. Using predicate = std::greater<>{} thus determines the first +// position in the sequence for which the elements is less than or equal to +// 'tolcrit'. Similarly, a predicate equivalent to '<=' returns the first +// position for which the elements is strictly greater than 'tolcrit'. + namespace { - struct RawTableEndPoints - { - struct { - std::vector gas; - std::vector water; - } connate; - - struct { - std::vector oil_in_gas; - std::vector oil_in_water; - std::vector gas; - std::vector water; - } critical; - - struct { - std::vector gas; - std::vector water; - } maximum; - }; + using ::Opm::satfunc::RawTableEndPoints; /* * See the "Saturation Functions" chapter in the Eclipse Technical @@ -116,8 +115,6 @@ namespace { return SatfuncFamily::none; } - // enum class limit { min, max }; - std::vector findMinWaterSaturation(const Opm::TableManager& tm, const Opm::Phases& ph) @@ -260,39 +257,67 @@ namespace { } } - /* - * These functions have been ported from an older implementation to instead - * use std::upper_bound and more from to make code -intent- - * clearer. This also made some (maybe intentional) details easier to spot. - * A short discussion: - * - * I don't know if not finding any element larger than 0.0 in the tables - * was ever supposed to happen (or even possible), in which case the vector - * elements remained at their initial value of 0.0. This behaviour has been - * preserved, but is now explicit. The original code was also not clear if - * it was possible to look up columns at index -1 (see critical_water for - * an example), but the new version is explicit about this. Unfortuately - * I'm not familiar enough with the maths or internal structure to make - * more than a guess here, but most of this behaviour should be preserved. - * - */ - - template - double critical_water( const T& table ) + template + auto crit_sat_index(const Opm::TableColumn& col, + const double tolcrit, + Predicate&& pred) { - const auto& col = table.getKrwColumn(); - const auto end = col.begin() + table.numRows(); - const auto critical = std::upper_bound( col.begin(), end, 0.0 ); - const auto index = std::distance( col.begin(), critical ); + using SizeT = std::remove_const_t< + std::remove_reference_t + >; - if( index == 0 || critical == end ) return 0.0; + auto begin = col.begin(); + auto pos = std::lower_bound(begin, col.end(), tolcrit, + std::forward(pred)); - return table.getSwColumn()[ index - 1 ]; + assert ((pos != col.end()) && + "Detected relative permeability function " + "without immobile state"); + + return static_cast(std::distance(begin, pos)); } - std::vector< double > + double crit_sat_increasing_KR(const Opm::TableColumn& sat, + const Opm::TableColumn& kr, + const double tolcrit) + { + // First position for which Kr(S) > tolcrit. + const auto i = crit_sat_index(kr, tolcrit, + [](const double kr1, const double kr2) + { + // kr1 <= kr2. Kr2 is 'tolcrit'. + return ! (kr2 < kr1); + }); + + return sat[i - 1]; // Last saturation for which Kr(S) <= tolcrit + } + + double crit_sat_decreasing_KR(const Opm::TableColumn& sat, + const Opm::TableColumn& kr, + const double tolcrit) + { + // First position for which Kr(S) <= tolcrit. + const auto i = crit_sat_index(kr, tolcrit, std::greater<>{}); + return sat[i]; + } + + /// Maximum water saturation for which Krw(Sw) <= tolcrit. + /// + /// Expected Table Format: + /// [Sw, Krw(Sw), ...other...] + /// + /// Krw increasing. + template + double critical_water(const T& table, const double tolcrit) + { + return crit_sat_increasing_KR(table.getSwColumn(), + table.getKrwColumn(), tolcrit); + } + + std::vector findCriticalWater(const Opm::TableManager& tm, - const Opm::Phases& ph) + const Opm::Phases& ph, + const double tolcrit) { const auto num_tables = tm.getTabdims().getNumSatTables(); @@ -302,12 +327,14 @@ namespace { const auto& swofTables = tm.getSwofTables(); const auto& swfnTables = tm.getSwfnTables(); - const auto famI = [&swofTables]( int i ) { - return critical_water( swofTables.getTable( i ) ); + const auto famI = [&swofTables, tolcrit](const int i) -> double + { + return critical_water(swofTables.getTable(i), tolcrit); }; - const auto famII = [&swfnTables]( int i ) { - return critical_water( swfnTables.getTable( i ) ); + const auto famII = [&swfnTables, tolcrit](const int i) -> double + { + return critical_water(swfnTables.getTable(i), tolcrit); }; switch( getSaturationFunctionFamily( tm, ph ) ) { @@ -317,31 +344,40 @@ namespace { } } - template< typename T > - double critical_gas( const T& table ) { - const auto& col = table.getKrgColumn(); - const auto end = col.begin() + table.numRows(); - const auto critical = std::upper_bound( col.begin(), end, 0.0 ); - const auto index = std::distance( col.begin(), critical ); - - if( index == 0 || critical == end ) return 0.0; - - return table.getSgColumn()[ index - 1 ]; + /// Maximum gas saturation for which Krg(Sg) <= tolcrit. + /// + /// Expected Table Format: + /// [Sg, Krg(Sg), ...other...] + /// + /// Krg increasing. + template + double critical_gas(const T& table, const double tolcrit) + { + return crit_sat_increasing_KR(table.getSgColumn(), + table.getKrgColumn(), tolcrit); } - double critical_gas( const Opm::SlgofTable& slgofTable ) { - const auto& col = slgofTable.getKrgColumn(); - const auto critical = std::upper_bound( col.begin(), col.end(), 0.0 ); - const auto index = std::distance( col.begin(), critical ); + /// Maximum gas saturation for which Krg(Sg) <= tolcrit. + /// + /// Table Format (Sl = So + Swco): + /// [Sl, Krg(Sl), Krog(Sl), Pcgo(Sl)] + /// + /// Krg decreasing, Krog increasing, Pcog not increasing. + double critical_gas(const Opm::SlgofTable& slgofTable, + const double tolcrit) + { + const auto sl_at_crit_gas = + crit_sat_decreasing_KR(slgofTable.getSlColumn(), + slgofTable.getKrgColumn(), tolcrit); - if( index == 0 || critical == col.end() ) return 0.0; - - return slgofTable.getSlColumn()[ index - 1 ]; + // Sg = 1 - Sl + return 1.0 - sl_at_crit_gas; } std::vector findCriticalGas(const Opm::TableManager& tm, - const Opm::Phases& ph) + const Opm::Phases& ph, + const double tolcrit) { const auto num_tables = tm.getTabdims().getNumSatTables(); @@ -352,16 +388,19 @@ namespace { const auto& sgofTables = tm.getSgofTables(); const auto& slgofTables = tm.getSlgofTables(); - const auto famI_sgof = [&sgofTables]( int i ) { - return critical_gas( sgofTables.getTable( i ) ); + const auto famI_sgof = [&sgofTables, tolcrit](const int i) -> double + { + return critical_gas(sgofTables.getTable(i), tolcrit); }; - const auto famI_slgof = [&slgofTables]( int i ) { - return critical_gas( slgofTables.getTable( i ) ); + const auto famI_slgof = [&slgofTables, tolcrit](const int i) -> double + { + return critical_gas(slgofTables.getTable(i), tolcrit); }; - const auto famII = [&sgfnTables]( int i ) { - return critical_gas( sgfnTables.getTable( i ) ); + const auto famII = [&sgfnTables, tolcrit](const int i) -> double + { + return critical_gas(sgfnTables.getTable(i), tolcrit); }; switch( getSaturationFunctionFamily( tm, ph ) ) { @@ -382,42 +421,53 @@ namespace { } } - double critical_oil_water( const Opm::SwofTable& swofTable ) { - const auto& col = swofTable.getKrowColumn(); + /// Maximum oil saturation for which Krow(So) <= tolcrit. + /// + /// Table Format: + /// [Sw, Krw(Sw), Krow(Sw), Pcow(Sw)] + /// + /// Krw increasing, Krow decreasing, Pcow not increasing. + double critical_oil_water(const Opm::SwofTable& swofTable, + const double tolcrit) + { + const auto sw_at_crit_oil = + crit_sat_decreasing_KR(swofTable.getSwColumn(), + swofTable.getKrowColumn(), tolcrit); - using reverse = std::reverse_iterator< decltype( col.begin() ) >; - auto rbegin = reverse( col.begin() + swofTable.numRows() ); - auto rend = reverse( col.begin() ); - const auto critical = std::upper_bound( rbegin, rend, 0.0 ); - const auto index = std::distance( col.begin(), critical.base() - 1 ); - - if( critical == rend ) return 0.0; - - return 1 - swofTable.getSwColumn()[ index + 1 ]; + // So = 1 - Sw + return 1.0 - sw_at_crit_oil; } - double critical_oil( const Opm::Sof2Table& sof2Table ) { - const auto& col = sof2Table.getKroColumn(); - const auto critical = std::upper_bound( col.begin(), col.end(), 0.0 ); - const auto index = std::distance( col.begin(), critical ); - - if( index == 0 || critical == col.end() ) return 0.0; - - return sof2Table.getSoColumn()[ index - 1 ]; + /// Maximum oil saturation for which Kro(So) <= tolcrit. + /// + /// Table Format: + /// [So, Kro(So)] + /// + /// Kro increasing. + double critical_oil(const Opm::Sof2Table& sof2Table, + const double tolcrit) + { + return crit_sat_increasing_KR(sof2Table.getSoColumn(), + sof2Table.getKroColumn(), tolcrit); } - double critical_oil( const Opm::Sof3Table& sof3Table, const Opm::TableColumn& col ) { - const auto critical = std::upper_bound( col.begin(), col.end(), 0.0 ); - const auto index = std::distance( col.begin(), critical ); - - if( index == 0 || critical == col.end() ) return 0.0; - - return sof3Table.getSoColumn()[ index - 1 ]; + /// Maximum oil saturation for which Kro(So) <= tolcrit. + /// + /// Table Format: + /// [So, Krow(So), Krog(So)] + /// + /// Krow increasing, Krog increasing. + double critical_oil(const Opm::Sof3Table& sof3Table, + const Opm::TableColumn& col, + const double tolcrit) + { + return crit_sat_increasing_KR(sof3Table.getSoColumn(), col, tolcrit); } std::vector findCriticalOilWater(const Opm::TableManager& tm, - const Opm::Phases& ph) + const Opm::Phases& ph, + const double tolcrit) { const auto num_tables = tm.getTabdims().getNumSatTables(); @@ -429,17 +479,20 @@ namespace { const auto& sof2Tables = tm.getSof2Tables(); const auto& sof3Tables = tm.getSof3Tables(); - const auto famI = [&swofTables]( int i ) { - return critical_oil_water( swofTables.getTable( i ) ); + const auto famI = [&swofTables, tolcrit](const int i) -> double + { + return critical_oil_water(swofTables.getTable(i), tolcrit); }; - const auto famII_2p = [&sof2Tables]( int i ) { - return critical_oil( sof2Tables.getTable( i ) ); + const auto famII_2p = [&sof2Tables, tolcrit](const int i) -> double + { + return critical_oil(sof2Tables.getTable(i), tolcrit); }; - const auto famII_3p = [&sof3Tables]( int i ) { - const auto& tb = sof3Tables.getTable( i ); - return critical_oil( tb, tb.getKrowColumn() ); + const auto famII_3p = [&sof3Tables, tolcrit](const int i) -> double + { + const auto& tb = sof3Tables.getTable(i); + return critical_oil(tb, tb.getKrowColumn(), tolcrit); }; switch( getSaturationFunctionFamily( tm, ph ) ) { @@ -453,36 +506,41 @@ namespace { } } - double critical_oil_gas( const Opm::SgofTable& sgofTable ) + /// Maximum oil saturation for which Krog(So) <= tolcrit. + /// + /// Table Format: + /// [Sg, Krg(Sg), Krog(Sg), Pcgo(Sg)] + /// + /// Krg increasing, Krog decreasing, Pcgo not decreasing. + double critical_oil_gas(const Opm::SgofTable& sgofTable, + const double tolcrit) { - const auto& col = sgofTable.getKrogColumn(); + const auto sg_at_crit_oil = + crit_sat_decreasing_KR(sgofTable.getSgColumn(), + sgofTable.getKrogColumn(), tolcrit); - using reverse = std::reverse_iterator< decltype( col.begin() ) >; - auto rbegin = reverse( col.begin() + sgofTable.numRows() ); - auto rend = reverse( col.begin() ); - const auto critical = std::upper_bound( rbegin, rend, 0.0 ); - if( critical == rend ) { - return 0.0; - } - const auto index = std::distance( col.begin(), critical.base() - 1 ); - return 1.0 - sgofTable.getSgColumn()[ index + 1 ]; + // So = 1 - Sg + return 1.0 - sg_at_crit_oil; } - double critical_oil_gas( const Opm::SlgofTable& sgofTable ) + /// Maximum oil saturation for which Krog(So) <= tolcrit. + /// + /// Table Format (Sl = So + Swco): + /// [Sl, Krg(Sl), Krog(Sl), Pcgo(Sl)] + /// + /// Krg decreasing, Krog increasing, Pcgo not increasing. + double critical_oil_gas(const Opm::SlgofTable& slgofTable, + const double tolcrit) { - const auto& col = sgofTable.getKrogColumn(); - const auto critical = std::upper_bound( col.begin(), col.end(), 0.0 ); - if (critical == col.end()) { - return 0.0; - } - const auto index = std::distance( col.begin(), critical - 1); - return sgofTable.getSlColumn()[ index ]; + return crit_sat_increasing_KR(slgofTable.getSlColumn(), + slgofTable.getKrogColumn(), tolcrit); } std::vector findCriticalOilGas(const Opm::TableManager& tm, const Opm::Phases& ph, - const std::vector& swco) + const std::vector& swco, + const double tolcrit) { const auto num_tables = tm.getTabdims().getNumSatTables(); @@ -495,23 +553,25 @@ namespace { const auto& sof2Tables = tm.getSof2Tables(); const auto& sof3Tables = tm.getSof3Tables(); - const auto famI_sgof = [&sgofTables, &swco](const int i) -> double + const auto famI_sgof = [&sgofTables, &swco, tolcrit](const int i) -> double { - return critical_oil_gas(sgofTables.getTable(i)) - swco[i]; + return critical_oil_gas(sgofTables.getTable(i), tolcrit) - swco[i]; }; - const auto famI_slgof = [&slgofTables, &swco](const int i) -> double + const auto famI_slgof = [&slgofTables, &swco, tolcrit](const int i) -> double { - return critical_oil_gas(slgofTables.getTable(i)) - swco[i]; + return critical_oil_gas(slgofTables.getTable(i), tolcrit) - swco[i]; }; - const auto famII_2p = [&sof2Tables]( int i ) { - return critical_oil( sof2Tables.getTable( i ) ); + const auto famII_2p = [&sof2Tables, tolcrit](const int i) -> double + { + return critical_oil(sof2Tables.getTable(i), tolcrit); }; - const auto famII_3p = [&sof3Tables]( int i ) { - const auto& tb = sof3Tables.getTable( i ); - return critical_oil( tb, tb.getKrogColumn() ); + const auto famII_3p = [&sof3Tables, tolcrit](const int i) -> double + { + const auto& tb = sof3Tables.getTable(i); + return critical_oil(tb, tb.getKrogColumn(), tolcrit); }; switch( getSaturationFunctionFamily( tm, ph ) ) { @@ -973,26 +1033,6 @@ namespace { } } - RawTableEndPoints - getRawTableEndpoints(const Opm::TableManager& tm, - const Opm::Phases& phases) - { - auto ep = RawTableEndPoints{}; - - ep.connate.gas = findMinGasSaturation(tm, phases); - ep.connate.water = findMinWaterSaturation(tm, phases); - - ep.critical.oil_in_gas = findCriticalOilGas(tm, phases, ep.connate.water); - ep.critical.oil_in_water = findCriticalOilWater(tm, phases); - ep.critical.gas = findCriticalGas(tm, phases); - ep.critical.water = findCriticalWater(tm, phases); - - ep.maximum.gas = findMaxGasSaturation(tm, phases); - ep.maximum.water = findMaxWaterSaturation(tm, phases); - - return ep; - } - double selectValue(const Opm::TableContainer& depthTables, int tableIdx, const std::string& columnName, @@ -1533,10 +1573,33 @@ namespace { } } // namespace Anonymous + +std::shared_ptr +Opm::satfunc::getRawTableEndpoints(const Opm::TableManager& tm, + const Opm::Phases& phases, + const double tolcrit) +{ + auto ep = std::make_shared(); + + ep->connate.gas = findMinGasSaturation(tm, phases); + ep->connate.water = findMinWaterSaturation(tm, phases); + + ep->critical.oil_in_gas = findCriticalOilGas(tm, phases, ep->connate.water, tolcrit); + ep->critical.oil_in_water = findCriticalOilWater(tm, phases, tolcrit); + ep->critical.gas = findCriticalGas(tm, phases, tolcrit); + ep->critical.water = findCriticalWater(tm, phases, tolcrit); + + ep->maximum.gas = findMaxGasSaturation(tm, phases); + ep->maximum.water = findMaxWaterSaturation(tm, phases); + + return ep; +} + std::vector Opm::satfunc::init(const std::string& keyword, const TableManager& tables, const Phases& phases, + const RawTableEndPoints& ep, const std::vector& cell_depth, const std::vector& num, const std::vector& endnum) @@ -1585,7 +1648,5 @@ Opm::satfunc::init(const std::string& keyword, + keyword + '\'' }; - const auto ep = getRawTableEndpoints(tables, phases); - return func->second(tables, phases, ep, cell_depth, num, endnum); } diff --git a/tests/parser/FieldPropsTests.cpp b/tests/parser/FieldPropsTests.cpp index 0065ac51b..7ac1771dd 100644 --- a/tests/parser/FieldPropsTests.cpp +++ b/tests/parser/FieldPropsTests.cpp @@ -19,8 +19,10 @@ #include #include +#include #include #include +#include #include #define BOOST_TEST_MODULE FieldPropsTests @@ -39,6 +41,7 @@ #include #include #include +#include #include #include "src/opm/parser/eclipse/EclipseState/Grid/FieldProps.hpp" @@ -588,6 +591,63 @@ MULTZ BOOST_CHECK_EQUAL( multx[0], 2 ); } +BOOST_AUTO_TEST_CASE(OPERATE) { + std::string deck_string = R"( +GRID + +PORO + 6*1.0 / + +OPERATE + PORO 1 3 1 1 1 1 'MAXLIM' PORO 0.50 / + PORO 1 3 2 2 1 1 'MAXLIM' PORO 0.25 / +/ + + +PERMX + 6*1/ + +PERMY + 6*1000/ + +OPERATE + PERMX 1 3 1 1 1 1 'MINLIM' PERMX 2 / + PERMX 1 3 2 2 1 1 'MINLIM' PERMX 4 / + PERMY 1 3 1 1 1 1 'MAXLIM' PERMY 100 / + PERMY 1 3 2 2 1 1 'MAXLIM' PERMY 200 / + PERMZ 1 3 1 1 1 1 'MULTA' PERMY 2 1000 / + PERMZ 1 3 2 2 1 1 'MULTA' PERMX 3 300 / +/ + + + + +)"; + + UnitSystem unit_system(UnitSystem::UnitType::UNIT_TYPE_METRIC); + auto to_si = [&unit_system](double raw_value) { return unit_system.to_si(UnitSystem::measure::permeability, raw_value); }; + EclipseGrid grid(3,2,1); + Deck deck = Parser{}.parseString(deck_string); + FieldPropsManager fpm(deck, Phases{true, true, true}, grid, TableManager()); + const auto& poro = fpm.get_double("PORO"); + BOOST_CHECK_EQUAL(poro[0], 0.50); + BOOST_CHECK_EQUAL(poro[3], 0.25); + + const auto& permx = fpm.get_double("PERMX"); + BOOST_CHECK_EQUAL(permx[0], to_si(2)); + BOOST_CHECK_EQUAL(permx[3], to_si(4)); + + const auto& permy = fpm.get_double("PERMY"); + BOOST_CHECK_EQUAL(permy[0], to_si(100)); + BOOST_CHECK_EQUAL(permy[3], to_si(200)); + + const auto& permz = fpm.get_double("PERMZ"); + for (std::size_t i = 0; i < 3; i++) { + BOOST_CHECK_EQUAL(permz[i] , 2*permy[i] + to_si(1000)); + BOOST_CHECK_EQUAL(permz[i+3], 3*permx[i+3] + to_si(300)); + } +} + namespace { std::string satfunc_model_setup() { @@ -869,6 +929,15 @@ SOF3 / )" }; } + std::string tolCrit(const double t) + { + std::ostringstream os; + + os << "TOLCRIT\n " << std::scientific << std::setw(11) << t << "\n/\n"; + + return os.str(); + } + std::string end() { return { R"( @@ -877,7 +946,370 @@ END } } // namespace Anonymous -BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I) { +BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_I_TolCrit_Zero) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_I() + end()) + }; + + const auto& tm = es.getTableManager(); + const auto& ph = es.runspec().phases(); + const auto tolcrit = 0.0; + + auto rtepPtr = satfunc::getRawTableEndpoints(tm, ph, tolcrit); + + // Water end-points + { + const auto swl = rtepPtr->connate .water; + const auto swcr = rtepPtr->critical.water; + const auto swu = rtepPtr->maximum .water; + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.071004, 1.0e-10); // == SWL. TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = rtepPtr->critical.oil_in_water; + const auto sogcr = rtepPtr->critical.oil_in_gas; + BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.791004, 1.0e-10); // TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(sogcr[0], 1.0 - 0.858996 - 0.071004, 1.0e-10); // Include SWL + } + + // Gas end-points + { + const auto sgl = rtepPtr->connate .gas; + const auto sgcr = rtepPtr->critical.gas; + const auto sgu = rtepPtr->maximum .gas; + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_II_TolCrit_Zero) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_II() + end()) + }; + + const auto& tm = es.getTableManager(); + const auto& ph = es.runspec().phases(); + const auto tolcrit = 0.0; + + auto rtepPtr = satfunc::getRawTableEndpoints(tm, ph, tolcrit); + + // Water end-points + { + const auto swl = rtepPtr->connate .water; + const auto swcr = rtepPtr->critical.water; + const auto swu = rtepPtr->maximum .water; + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.071004, 1.0e-10); // == SWL. TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = rtepPtr->critical.oil_in_water; + const auto sogcr = rtepPtr->critical.oil_in_gas; + BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.791004, 1.0e-10); // TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(sogcr[0], 1.0 - 0.858996 - 0.071004, 1.0e-10); // Include SWL + } + + // Gas end-points + { + const auto sgl = rtepPtr->connate .gas; + const auto sgcr = rtepPtr->critical.gas; + const auto sgu = rtepPtr->maximum .gas; + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_I_TolCrit_Default) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_I() + end()) + }; + + const auto& tm = es.getTableManager(); + const auto& ph = es.runspec().phases(); + const auto tolcrit = es.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); // 1.0e-6. + + auto rtepPtr = satfunc::getRawTableEndpoints(tm, ph, tolcrit); + + // Water end-points + { + const auto swl = rtepPtr->connate .water; + const auto swcr = rtepPtr->critical.water; + const auto swu = rtepPtr->maximum .water; + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.091004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = rtepPtr->critical.oil_in_water; + const auto sogcr = rtepPtr->critical.oil_in_gas; + BOOST_CHECK_CLOSE(sowcr[0], 0.228996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.070000, 1.0e-10); // Max So for which Krog(So) <= TOLCRIT + } + + // Gas end-points + { + const auto sgl = rtepPtr->connate .gas; + const auto sgcr = rtepPtr->critical.gas; + const auto sgu = rtepPtr->maximum .gas; + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); // Max Sg for which Krg(Sg) <= TOLCRIT + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_II_TolCrit_Default) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_II() + end()) + }; + + const auto& tm = es.getTableManager(); + const auto& ph = es.runspec().phases(); + const auto tolcrit = es.runspec().saturationFunctionControls() + .minimumRelpermMobilityThreshold(); // 1.0e-6. + + auto rtepPtr = satfunc::getRawTableEndpoints(tm, ph, tolcrit); + + // Water end-points + { + const auto swl = rtepPtr->connate .water; + const auto swcr = rtepPtr->critical.water; + const auto swu = rtepPtr->maximum .water; + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.091004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = rtepPtr->critical.oil_in_water; + const auto sogcr = rtepPtr->critical.oil_in_gas; + BOOST_CHECK_CLOSE(sowcr[0], 0.228996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.070000, 1.0e-10); // Max So for which Krog(So) <= TOLCRIT + } + + // Gas end-points + { + const auto sgl = rtepPtr->connate .gas; + const auto sgcr = rtepPtr->critical.gas; + const auto sgu = rtepPtr->maximum .gas; + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); // Max Sg for which Krg(Sg) <= TOLCRIT + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_I_TolCrit_Large) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_I() + end()) + }; + + const auto& tm = es.getTableManager(); + const auto& ph = es.runspec().phases(); + const auto tolcrit = 0.01; + + auto rtepPtr = satfunc::getRawTableEndpoints(tm, ph, tolcrit); + + // Water end-points + { + const auto swl = rtepPtr->connate .water; + const auto swcr = rtepPtr->critical.water; + const auto swu = rtepPtr->maximum .water; + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.231004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = rtepPtr->critical.oil_in_water; + const auto sogcr = rtepPtr->critical.oil_in_gas; + BOOST_CHECK_CLOSE(sowcr[0], 0.448996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.238996, 1.0e-10); // Max So for which Krog(So) <= TOLCRIT + } + + // Gas end-points + { + const auto sgl = rtepPtr->connate .gas; + const auto sgcr = rtepPtr->critical.gas; + const auto sgu = rtepPtr->maximum .gas; + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.09, 1.0e-10); // Max Sg for which Krg(Sg) <= TOLCRIT + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_II_TolCrit_Large) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_II() + end()) + }; + + const auto& tm = es.getTableManager(); + const auto& ph = es.runspec().phases(); + const auto tolcrit = 0.01; + + auto rtepPtr = satfunc::getRawTableEndpoints(tm, ph, tolcrit); + + // Water end-points + { + const auto swl = rtepPtr->connate .water; + const auto swcr = rtepPtr->critical.water; + const auto swu = rtepPtr->maximum .water; + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.231004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = rtepPtr->critical.oil_in_water; + const auto sogcr = rtepPtr->critical.oil_in_gas; + BOOST_CHECK_CLOSE(sowcr[0], 0.448996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.248996, 1.0e-10); // Max So for which Krog(So) <= TOLCRIT + } + + // Gas end-points + { + const auto sgl = rtepPtr->connate .gas; + const auto sgcr = rtepPtr->critical.gas; + const auto sgu = rtepPtr->maximum .gas; + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.09, 1.0e-10); // Max Sg for which Krg(Sg) <= TOLCRIT + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + } +} + +// ===================================================================== + +BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I_TolCrit_Zero) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + tolCrit(0.0) + satfunc_family_I() + end()) + }; + + auto fp = es.fieldProps(); + + // Water end-points + { + const auto swl = fp.get_double("SWL"); + const auto swcr = fp.get_double("SWCR"); + const auto swu = fp.get_double("SWU"); + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.071004, 1.0e-10); // == SWL. TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + + const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.791004) + const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1) + BOOST_CHECK_CLOSE(krwr[0], 0.911429, 1.0e-10); + BOOST_CHECK_CLOSE(krw [0], 1.0 , 1.0e-10); + + const auto pcw = fp.get_double("PCW"); + BOOST_CHECK_CLOSE(pcw[0], 7.847999*unit::barsa, 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = fp.get_double("SOWCR"); + const auto sogcr = fp.get_double("SOGCR"); + BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.791004, 1.0e-10); // TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(sogcr[0], 1.0 - 0.858996 - 0.071004, 1.0e-10); // Include SWL + + const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.928996) + const auto krorg = fp.get_double("KRORG"); // Krog(So=1-Sgcr-Swl) = Krog(So=0.898996) + const auto kro = fp.get_double("KRO"); // Krow(So=Somax) = Krog(So=Somax) + BOOST_CHECK_CLOSE(krorw[0], 1.0, 1.0e-10); // TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(krorg[0], 0.896942, 1.0e-10); + BOOST_CHECK_CLOSE(kro [0], 1.0, 1.0e-10); + } + + // Gas end-points + { + const auto sgl = fp.get_double("SGL"); + const auto sgcr = fp.get_double("SGCR"); + const auto sgu = fp.get_double("SGU"); + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + + const auto krgr = fp.get_double("KRGR"); // Krg(Sg=1-Sogcr-Swl) = Krg(Sg=0.858996) + const auto krg = fp.get_double("KRG"); // Krg(Sgmax) = Krg(Sg=0.928996) + BOOST_CHECK_CLOSE(krgr[0], 0.866135, 1.0e-10); + BOOST_CHECK_CLOSE(krg [0], 1.0 , 1.0e-10); + + const auto pcg = fp.get_double("PCG"); + BOOST_CHECK_CLOSE(pcg[0], 0.0, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II_TolCrit_Zero) { + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + tolCrit(0.0) + satfunc_family_II() + end()) + }; + + auto fp = es.fieldProps(); + + // Water end-points + { + const auto swl = fp.get_double("SWL"); + const auto swcr = fp.get_double("SWCR"); + const auto swu = fp.get_double("SWU"); + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.071004, 1.0e-10); // == SWL. TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + + const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.791004) + const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1) + BOOST_CHECK_CLOSE(krwr[0], 0.911429, 1.0e-10); + BOOST_CHECK_CLOSE(krw [0], 1.0 , 1.0e-10); + + const auto pcw = fp.get_double("PCW"); + BOOST_CHECK_CLOSE(pcw[0], 7.847999*unit::barsa, 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = fp.get_double("SOWCR"); + const auto sogcr = fp.get_double("SOGCR"); + BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.791004, 1.0e-10); // TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(sogcr[0], 1.0 - 0.858996 - 0.071004, 1.0e-10); // Include SWL + + const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.928996) + const auto krorg = fp.get_double("KRORG"); // Krog(So=1-Sgcr-Swl) = Krog(So=0.898996) + const auto kro = fp.get_double("KRO"); // Krow(So=Somax) = Krog(So=Somax) + BOOST_CHECK_CLOSE(krorw[0], 1.0, 1.0e-10); // TOLCRIT = 0.0 + BOOST_CHECK_CLOSE(krorg[0], 0.896942, 1.0e-10); + BOOST_CHECK_CLOSE(kro [0], 1.0, 1.0e-10); + } + + // Gas end-points + { + const auto sgl = fp.get_double("SGL"); + const auto sgcr = fp.get_double("SGCR"); + const auto sgu = fp.get_double("SGU"); + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + + const auto krgr = fp.get_double("KRGR"); // Krg(Sg=1-Sogcr-Swl) = Krg(Sg=0.858996) + const auto krg = fp.get_double("KRG"); // Krg(Sgmax) = Krg(Sg=0.928996) + BOOST_CHECK_CLOSE(krgr[0], 0.866135, 1.0e-10); + BOOST_CHECK_CLOSE(krg [0], 1.0 , 1.0e-10); + + const auto pcg = fp.get_double("PCG"); + BOOST_CHECK_CLOSE(pcg[0], 0.0, 1.0e-10); + } +} + +BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I_TolCrit_Default) { + // TOLCRIT = 1.0e-6 const auto es = ::Opm::EclipseState { ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_I() + end()) }; @@ -890,12 +1322,12 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I) { const auto swcr = fp.get_double("SWCR"); const auto swu = fp.get_double("SWU"); BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); - BOOST_CHECK_CLOSE(swcr[0], 0.071004, 1.0e-10); // == SWL. TOLCRIT not currently implemented. + BOOST_CHECK_CLOSE(swcr[0], 0.091004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); - const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.791004) + const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.771004) const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1) - BOOST_CHECK_CLOSE(krwr[0], 0.911429, 1.0e-10); + BOOST_CHECK_CLOSE(krwr[0], 0.835916, 1.0e-10); BOOST_CHECK_CLOSE(krw [0], 1.0 , 1.0e-10); const auto pcw = fp.get_double("PCW"); @@ -906,13 +1338,13 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I) { { const auto sowcr = fp.get_double("SOWCR"); const auto sogcr = fp.get_double("SOGCR"); - BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.791004, 1.0e-10); // TOLCRIT currently not supported. + BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.771004, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT. BOOST_CHECK_CLOSE(sogcr[0], 1.0 - 0.858996 - 0.071004, 1.0e-10); // Include SWL - const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.928996) + const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.908996) const auto krorg = fp.get_double("KRORG"); // Krog(So=1-Sgcr-Swl) = Krog(So=0.898996) const auto kro = fp.get_double("KRO"); // Krow(So=Somax) = Krog(So=Somax) - BOOST_CHECK_CLOSE(krorw[0], 1.0, 1.0e-10); // TOLCRIT currently not supported. + BOOST_CHECK_CLOSE(krorw[0], 0.882459, 1.0e-10); BOOST_CHECK_CLOSE(krorg[0], 0.896942, 1.0e-10); BOOST_CHECK_CLOSE(kro [0], 1.0, 1.0e-10); } @@ -936,7 +1368,8 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I) { } } -BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II) { +BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II_TolCrit_Default) { + // TOLCRIT = 1.0e-6 const auto es = ::Opm::EclipseState { ::Opm::Parser{}.parseString(satfunc_model_setup() + satfunc_family_II() + end()) }; @@ -949,12 +1382,12 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II) { const auto swcr = fp.get_double("SWCR"); const auto swu = fp.get_double("SWU"); BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); - BOOST_CHECK_CLOSE(swcr[0], 0.071004, 1.0e-10); // == SWL. TOLCRIT not currently implemented. + BOOST_CHECK_CLOSE(swcr[0], 0.091004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); - const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.791004) + const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.771004) const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1) - BOOST_CHECK_CLOSE(krwr[0], 0.911429, 1.0e-10); + BOOST_CHECK_CLOSE(krwr[0], 0.835916, 1.0e-10); BOOST_CHECK_CLOSE(krw [0], 1.0 , 1.0e-10); const auto pcw = fp.get_double("PCW"); @@ -965,21 +1398,21 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II) { { const auto sowcr = fp.get_double("SOWCR"); const auto sogcr = fp.get_double("SOGCR"); - BOOST_CHECK_CLOSE(sowcr[0], 1.0 - 0.791004, 1.0e-10); // TOLCRIT currently not supported. - BOOST_CHECK_CLOSE(sogcr[0], 1.0 - 0.858996 - 0.071004, 1.0e-10); // Include SWL + BOOST_CHECK_CLOSE(sowcr[0], 0.228996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.070000, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT - const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.928996) + const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.908996) const auto krorg = fp.get_double("KRORG"); // Krog(So=1-Sgcr-Swl) = Krog(So=0.898996) const auto kro = fp.get_double("KRO"); // Krow(So=Somax) = Krog(So=Somax) - BOOST_CHECK_CLOSE(krorw[0], 1.0, 1.0e-10); // TOLCRIT currently not supported. + BOOST_CHECK_CLOSE(krorw[0], 0.882459, 1.0e-10); BOOST_CHECK_CLOSE(krorg[0], 0.896942, 1.0e-10); - BOOST_CHECK_CLOSE(kro [0], 1.0, 1.0e-10); + BOOST_CHECK_CLOSE(kro [0], 1.0 , 1.0e-10); } // Gas end-points { const auto sgl = fp.get_double("SGL"); - const auto sgcr = fp.get_double("SGCR"); + const auto sgcr = fp.get_double("SGCR"); // Max Sg for which Krg(Sg) <= TOLCRIT const auto sgu = fp.get_double("SGU"); BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); BOOST_CHECK_CLOSE(sgcr[0], 0.03, 1.0e-10); @@ -995,61 +1428,122 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II) { } } -BOOST_AUTO_TEST_CASE(OPERATE) { - std::string deck_string = R"( -GRID +BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_I_TolCrit_Large) { + // TOLCRIT = 0.01 + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + tolCrit(0.01) + satfunc_family_I() + end()) + }; -PORO - 6*1.0 / + auto fp = es.fieldProps(); -OPERATE - PORO 1 3 1 1 1 1 'MAXLIM' PORO 0.50 / - PORO 1 3 2 2 1 1 'MAXLIM' PORO 0.25 / -/ + // Water end-points + { + const auto swl = fp.get_double("SWL"); + const auto swcr = fp.get_double("SWCR"); + const auto swu = fp.get_double("SWU"); + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.231004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.551004) + const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1) + BOOST_CHECK_CLOSE(krwr[0], 0.261115, 1.0e-10); + BOOST_CHECK_CLOSE(krw [0], 1.0 , 1.0e-10); -PERMX - 6*1/ + const auto pcw = fp.get_double("PCW"); + BOOST_CHECK_CLOSE(pcw[0], 7.847999*unit::barsa, 1.0e-10); + } -PERMY - 6*1000/ + // Oil end-points + { + const auto sowcr = fp.get_double("SOWCR"); + const auto sogcr = fp.get_double("SOGCR"); + BOOST_CHECK_CLOSE(sowcr[0], 0.448996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.238996, 1.0e-10); // Max So for which Krog(So) <= TOLCRIT -OPERATE - PERMX 1 3 1 1 1 1 'MINLIM' PERMX 2 / - PERMX 1 3 2 2 1 1 'MINLIM' PERMX 4 / - PERMY 1 3 1 1 1 1 'MAXLIM' PERMY 100 / - PERMY 1 3 2 2 1 1 'MAXLIM' PERMY 200 / - PERMZ 1 3 1 1 1 1 'MULTA' PERMY 2 1000 / - PERMZ 1 3 2 2 1 1 'MULTA' PERMX 3 300 / -/ + const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.768996) + const auto krorg = fp.get_double("KRORG"); // Krog(So=1-Sgcr-Swl) = Krog(So=0.838996) + const auto kro = fp.get_double("KRO"); // Krow(So=Somax) = Krog(So=Somax) + BOOST_CHECK_CLOSE(krorw[0], 0.328347, 1.0e-10); + BOOST_CHECK_CLOSE(krorg[0], 0.712749, 1.0e-10); + BOOST_CHECK_CLOSE(kro [0], 1.0, 1.0e-10); + } + // Gas end-points + { + const auto sgl = fp.get_double("SGL"); + const auto sgcr = fp.get_double("SGCR"); + const auto sgu = fp.get_double("SGU"); + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.090000, 1.0e-10); // Max Sg for which Krg(Sg) <= TOLCRIT + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + const auto krgr = fp.get_double("KRGR"); // Krg(Sg=1-Sogcr-Swl) = Krg(Sg=0.690000) + const auto krg = fp.get_double("KRG"); // Krg(Sgmax) = Krg(Sg=0.928996) + BOOST_CHECK_CLOSE(krgr[0], 0.578171, 1.0e-10); + BOOST_CHECK_CLOSE(krg [0], 1.0 , 1.0e-10); - -)"; - - UnitSystem unit_system(UnitSystem::UnitType::UNIT_TYPE_METRIC); - auto to_si = [&unit_system](double raw_value) { return unit_system.to_si(UnitSystem::measure::permeability, raw_value); }; - EclipseGrid grid(3,2,1); - Deck deck = Parser{}.parseString(deck_string); - FieldPropsManager fpm(deck, Phases{true, true, true}, grid, TableManager()); - const auto& poro = fpm.get_double("PORO"); - BOOST_CHECK_EQUAL(poro[0], 0.50); - BOOST_CHECK_EQUAL(poro[3], 0.25); - - const auto& permx = fpm.get_double("PERMX"); - BOOST_CHECK_EQUAL(permx[0], to_si(2)); - BOOST_CHECK_EQUAL(permx[3], to_si(4)); - - const auto& permy = fpm.get_double("PERMY"); - BOOST_CHECK_EQUAL(permy[0], to_si(100)); - BOOST_CHECK_EQUAL(permy[3], to_si(200)); - - const auto& permz = fpm.get_double("PERMZ"); - for (std::size_t i = 0; i < 3; i++) { - BOOST_CHECK_EQUAL(permz[i] , 2*permy[i] + to_si(1000)); - BOOST_CHECK_EQUAL(permz[i+3], 3*permx[i+3] + to_si(300)); + const auto pcg = fp.get_double("PCG"); + BOOST_CHECK_CLOSE(pcg[0], 0.0, 1.0e-10); } } +BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II_TolCrit_Large) { + // TOLCRIT = 0.01 + const auto es = ::Opm::EclipseState { + ::Opm::Parser{}.parseString(satfunc_model_setup() + tolCrit(0.01) + satfunc_family_II() + end()) + }; + auto fp = es.fieldProps(); + + // Water end-points + { + const auto swl = fp.get_double("SWL"); + const auto swcr = fp.get_double("SWCR"); + const auto swu = fp.get_double("SWU"); + BOOST_CHECK_CLOSE(swl [0], 0.071004, 1.0e-10); + BOOST_CHECK_CLOSE(swcr[0], 0.231004, 1.0e-10); // Max Sw for which Krw(Sw) <= TOLCRIT + BOOST_CHECK_CLOSE(swu [0], 1.0 , 1.0e-10); + + const auto krwr = fp.get_double("KRWR"); // Krw(Sw=1-Sowcr-Sgl) = Krw(Sw=0.551004) + const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1) + BOOST_CHECK_CLOSE(krwr[0], 0.261115, 1.0e-10); + BOOST_CHECK_CLOSE(krw [0], 1.0 , 1.0e-10); + + const auto pcw = fp.get_double("PCW"); + BOOST_CHECK_CLOSE(pcw[0], 7.847999*unit::barsa, 1.0e-10); + } + + // Oil end-points + { + const auto sowcr = fp.get_double("SOWCR"); + const auto sogcr = fp.get_double("SOGCR"); + BOOST_CHECK_CLOSE(sowcr[0], 0.448996, 1.0e-10); // Max So for which Krow(So) <= TOLCRIT + BOOST_CHECK_CLOSE(sogcr[0], 0.248996, 1.0e-10); // Max So for which Krog(So) <= TOLCRIT + + const auto krorw = fp.get_double("KRORW"); // Krow(So=1-Swcr-Sgl) = Krow(So=0.768996) + const auto krorg = fp.get_double("KRORG"); // Krog(So=1-Sgcr-Swl) = Krog(So=0.838996) + const auto kro = fp.get_double("KRO"); // Krow(So=Somax) = Krog(So=Somax) + BOOST_CHECK_CLOSE(krorw[0], 0.328347, 1.0e-10); + BOOST_CHECK_CLOSE(krorg[0], 0.712749, 1.0e-10); + BOOST_CHECK_CLOSE(kro [0], 1.0, 1.0e-10); + } + + // Gas end-points + { + const auto sgl = fp.get_double("SGL"); + const auto sgcr = fp.get_double("SGCR"); + const auto sgu = fp.get_double("SGU"); + BOOST_CHECK_CLOSE(sgl [0], 0.0, 1.0e-10); + BOOST_CHECK_CLOSE(sgcr[0], 0.090000, 1.0e-10); // Max Sg for which Krg(Sg) <= TOLCRIT + BOOST_CHECK_CLOSE(sgu [0], 0.928996, 1.0e-10); + + const auto krgr = fp.get_double("KRGR"); // Krg(Sg=1-Sogcr-Swl) = Krg(Sg=0.680000) + const auto krg = fp.get_double("KRG"); // Krg(Sgmax) = Krg(Sg=0.928996) + BOOST_CHECK_CLOSE(krgr[0], 0.562914, 1.0e-10); + BOOST_CHECK_CLOSE(krg [0], 1.0 , 1.0e-10); + + const auto pcg = fp.get_double("PCG"); + BOOST_CHECK_CLOSE(pcg[0], 0.0, 1.0e-10); + } +}