OilfieldToolsNet/opm-common
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Sat-Func Initializers: Add TOLCRIT Support

Browse Source 
This commit makes the 'critical_*()' helper functions aware of the
TOLCRIT feature.  We now consider a phase to be immobile at a one
saturation value if its relative permeability at that saturation is
less than or equal to TOLCRIT.  TOLCRIT mirrors E100's default value
of 1.0e-6 unless otherwise specified in the PROPS section.

All table scanners are implemented in terms of the expression

    std::lower_bound(begin, end, tolcrit, predicate)

which ultimately 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'.

Add a set of unit tests to demonstrate and exercise the new feature.
This commit is contained in:
Bård Skaflestad
2020-06-10 14:33:01 +02:00
parent bcdc8bb119
commit abec96a387
2 changed files with 663 additions and 110 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/opm/parser/eclipse/EclipseState/Grid/SatfuncPropertyInitializers.cpp
+176 -65
View File
@@ -28,13 +28,16 @@
#include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/Tabdims.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableContainer.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableColumn.hpp>
#include <opm/parser/eclipse/EclipseState/Tables/TableManager.hpp>
#include <opm/parser/eclipse/Utility/Functional.hpp>
#include <algorithm>
#include <array>
#include <cassert>
#include <exception>
#include <functional>
#include <iterator>
#include <memory>
#include <stdexcept>
@@ -42,6 +45,17 @@
#include <stddef.h>
// 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 {
using ::Opm::satfunc::RawTableEndPoints;
@@ -243,37 +257,36 @@ namespace {
}
}
/*
* These functions have been ported from an older implementation to instead
* use std::upper_bound and more from <algorithm> 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.
*
*/
/// Maximum water saturation for which Krw(Sw) <= tolcrit.
///
/// Expected Table Format:
/// [Sw, Krw(Sw), ...other...]
///
/// Krw increasing.
template <typename T>
double critical_water(const T& table, const double tolcrit)
{
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 );
auto begin = col.begin();
auto end = col.end();
if( index == 0 || critical == end ) return 0.0;
// Locate first position for which Krw(Sw) > tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit,
[](const double kr1, const double kr2) -> bool
{
// kr1 <= kr2. Kr2 is 'tolcrit'.
return ! (kr2 < kr1);
});
return table.getSwColumn()[ index - 1 ];
assert ((critPos != begin) &&
"Detected Krw(Sw) Without Immobile Water Saturation");
const auto index = std::distance(begin, critPos);
return table.getSwColumn()[index - 1];
}
std::vector< double >
std::vector<double>
findCriticalWater(const Opm::TableManager& tm,
const Opm::Phases& ph,
const double tolcrit)
@@ -303,29 +316,58 @@ namespace {
}
}
template< typename T >
/// Maximum gas saturation for which Krg(Sg) <= tolcrit.
///
/// Expected Table Format:
/// [Sg, Krg(Sg), ...other...]
///
/// Krg increasing.
template <typename T>
double critical_gas(const T& table, const double tolcrit)
{
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 );
auto begin = col.begin();
auto end = col.end();
if( index == 0 || critical == end ) return 0.0;
// Locate first position for which Krg(Sg) > tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit,
[](const double kr1, const double kr2) -> bool
{
// kr1 <= kr2. Kr2 is 'tolcrit'.
return ! (kr2 < kr1);
});
return table.getSgColumn()[ index - 1 ];
assert ((critPos != begin) &&
"Detected Krg(Sg) Without Immobile Gas Saturation");
const auto index = std::distance(begin, critPos);
return table.getSgColumn()[index - 1];
}
/// 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& col = slgofTable.getKrgColumn();
const auto critical = std::upper_bound( col.begin(), col.end(), 0.0 );
const auto index = std::distance( col.begin(), critical );
auto begin = col.begin();
auto end = col.end();
if( index == 0 || critical == col.end() ) return 0.0;
// Locate first position for which Krg(Sl) <= tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit, std::greater<>{});
return slgofTable.getSlColumn()[ index - 1 ];
assert ((critPos != end) &&
"Detected Krg(Sl) Without Immobile Gas Saturation");
const auto index = std::distance(begin, critPos);
// Sg = 1 - Sl
return 1.0 - slgofTable.getSlColumn()[index]; // Not "index - 1"
}
std::vector<double>
@@ -375,44 +417,87 @@ namespace {
}
}
/// 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& col = swofTable.getKrowColumn();
auto begin = col.begin();
auto end = col.end();
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 );
// Locate first position for which Krow(So) <= tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit, std::greater<>{});
if( critical == rend ) return 0.0;
assert ((critPos != end) &&
"Detected Krow(So) Without Immobile Oil Saturation");
return 1 - swofTable.getSwColumn()[ index + 1 ];
const auto index = std::distance(begin, critPos);
// So = 1 - Sw
return 1.0 - swofTable.getSwColumn()[index]; // Not "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)
{
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 );
auto begin = col.begin();
auto end = col.end();
if( index == 0 || critical == col.end() ) return 0.0;
// Locate first position for which Kro(So) > tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit,
[](const double kr1, const double kr2) -> bool
{
// kr1 <= kr2. Kr2 is 'tolcrit'.
return ! (kr2 < kr1);
});
return sof2Table.getSoColumn()[ index - 1 ];
assert ((critPos != begin) &&
"Detected Kro(So) Without Immobile Oil Saturation");
const auto index = std::distance(begin, critPos);
return sof2Table.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)
{
const auto critical = std::upper_bound( col.begin(), col.end(), 0.0 );
const auto index = std::distance( col.begin(), critical );
auto begin = col.begin();
auto end = col.end();
if( index == 0 || critical == col.end() ) return 0.0;
// Locate first position for which Kro(So) > tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit,
[](const double kr1, const double kr2) -> bool
{
// kr1 <= kr2. Kr2 is 'tolcrit'.
return ! (kr2 < kr1);
});
return sof3Table.getSoColumn()[ index - 1 ];
assert ((critPos != begin) &&
"Detected Kro(So) Without Immobile Oil Saturation");
const auto index = std::distance(begin, critPos);
return sof3Table.getSoColumn()[index - 1];
}
std::vector<double>
@@ -457,32 +542,58 @@ namespace {
}
}
/// 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();
auto begin = col.begin();
auto end = col.end();
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 ];
// Locate first position for which Krog(So) <= tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit, std::greater<>{});
assert ((critPos != end) &&
"Detected Krog(So) Without Immobile Oil Saturation");
const auto index = std::distance(begin, critPos);
// So = 1 - Sg
return 1.0 - sgofTable.getSgColumn()[index]; // Not "index - 1"
}
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 ];
const auto& col = slgofTable.getKrogColumn();
auto begin = col.begin();
auto end = col.end();
// Locate first position for which Krog(So) > tolcrit.
auto critPos = std::lower_bound(begin, end, tolcrit,
[](const double kr1, const double kr2) -> bool
{
// kr1 <= kr2. Kr2 is tolcrit'.
return ! (kr2 < kr1);
});
assert ((critPos != end) &&
"Detected Krog(So) Without Immobile Oil Saturation");
const auto index = std::distance(begin, critPos);
return slgofTable.getSlColumn()[index - 1];
}
std::vector<double>
tests/parser/FieldPropsTests.cpp
+487 -45
View File
@@ -591,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()
{
@@ -929,6 +986,210 @@ BOOST_AUTO_TEST_CASE(RawTableEndPoints_Family_I_TolCrit_Zero) {
}
}
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())
@@ -1047,61 +1308,242 @@ BOOST_AUTO_TEST_CASE(SatFunc_EndPts_Family_II_TolCrit_Zero) {
}
}
BOOST_AUTO_TEST_CASE(OPERATE) {
std::string deck_string = R"(
GRID
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())
};
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.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.771004)
const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1)
BOOST_CHECK_CLOSE(krwr[0], 0.835916, 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], 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
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.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], 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);
}
// 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);
)";
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_Default) {
// TOLCRIT = 1.0e-6
const auto es = ::Opm::EclipseState {
::Opm::Parser{}.parseString(satfunc_model_setup() + 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.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.771004)
const auto krw = fp.get_double("KRW"); // Krw(Swmax) = Krw(Sw=1)
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");
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.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.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], 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);
}
// Gas end-points
{
const auto sgl = fp.get_double("SGL");
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);
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_Large) {
// TOLCRIT = 0.01
const auto es = ::Opm::EclipseState {
::Opm::Parser{}.parseString(satfunc_model_setup() + tolCrit(0.01) + 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.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.238996, 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.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);
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);
}
}