diff --git a/opm/autodiff/VFPProperties.hpp b/opm/autodiff/VFPProperties.hpp index 11f484852..f6c468903 100644 --- a/opm/autodiff/VFPProperties.hpp +++ b/opm/autodiff/VFPProperties.hpp @@ -70,29 +70,29 @@ namespace Opm { VFPProperties(int table_num, double datum_depth, - FLO_TYPE flo_type, - WFR_TYPE wfr_type, - GFR_TYPE gfr_type, - ALQ_TYPE alq_type, - const std::vector& flo_data, - const std::vector& thp_data, - const std::vector& wfr_data, - const std::vector& gfr_data, - const std::vector& alq_data, - array_type data - ) : - table_num_(table_num), - datum_depth_(datum_depth), - flo_type_(flo_type), - wfr_type_(wfr_type), - gfr_type_(gfr_type), - alq_type_(alq_type), - flo_data_(flo_data), - thp_data_(thp_data), - wfr_data_(wfr_data), - gfr_data_(gfr_data), - alq_data_(alq_data), - data_(data) { + FLO_TYPE flo_type, + WFR_TYPE wfr_type, + GFR_TYPE gfr_type, + ALQ_TYPE alq_type, + const std::vector& flo_data, + const std::vector& thp_data, + const std::vector& wfr_data, + const std::vector& gfr_data, + const std::vector& alq_data, + array_type data + ) : + table_num_(table_num), + datum_depth_(datum_depth), + flo_type_(flo_type), + wfr_type_(wfr_type), + gfr_type_(gfr_type), + alq_type_(alq_type), + flo_data_(flo_data), + thp_data_(thp_data), + wfr_data_(wfr_data), + gfr_data_(gfr_data), + alq_data_(alq_data), + data_(data) { } @@ -105,99 +105,99 @@ namespace Opm { //Rate type try { - std::string flo_string = header->getItem("RATE_TYPE")->getString(0); - if (flo_string == "OIL") { - flo_type_ = FLO_OIL; - } - else if (flo_string == "LIQ") { - flo_type_ = FLO_LIQ; - } - else if (flo_string == "GAS") { - flo_type_ = FLO_GAS; - } - else { - flo_type_ = FLO_INVALID; - } + std::string flo_string = header->getItem("RATE_TYPE")->getString(0); + if (flo_string == "OIL") { + flo_type_ = FLO_OIL; + } + else if (flo_string == "LIQ") { + flo_type_ = FLO_LIQ; + } + else if (flo_string == "GAS") { + flo_type_ = FLO_GAS; + } + else { + flo_type_ = FLO_INVALID; + } + } + catch (std::invalid_argument& e) { + //TODO: log here + flo_type_ = FLO_INVALID; } - catch (std::invalid_argument& e) { - //TODO: log here - flo_type_ = FLO_INVALID; - } //Water fraction try { - std::string wfr_string = header->getItem("WFR")->getString(0); - if (wfr_string == "WOR") { - wfr_type_ = WFR_WOR; - } - else if (wfr_string == "WCT") { - wfr_type_ = WFR_WCT; - } - else if (wfr_string == "WGR") { - wfr_type_ = WFR_WGR; - } - else { - wfr_type_ = WFR_INVALID; - } + std::string wfr_string = header->getItem("WFR")->getString(0); + if (wfr_string == "WOR") { + wfr_type_ = WFR_WOR; + } + else if (wfr_string == "WCT") { + wfr_type_ = WFR_WCT; + } + else if (wfr_string == "WGR") { + wfr_type_ = WFR_WGR; + } + else { + wfr_type_ = WFR_INVALID; + } + } + catch (std::invalid_argument& e) { + //TODO: log here + wfr_type_ = WFR_INVALID; } - catch (std::invalid_argument& e) { - //TODO: log here - wfr_type_ = WFR_INVALID; - } //Gas fraction try { - std::string gfr_string = header->getItem("GFR")->getString(0); - if (gfr_string == "GOR") { - gfr_type_ = GFR_GOR; - } - else if (gfr_string == "GLR") { - gfr_type_ = GFR_GLR; - } - else if (gfr_string == "OGR") { - gfr_type_ = GFR_OGR; - } - else { - gfr_type_ = GFR_INVALID; - } + std::string gfr_string = header->getItem("GFR")->getString(0); + if (gfr_string == "GOR") { + gfr_type_ = GFR_GOR; + } + else if (gfr_string == "GLR") { + gfr_type_ = GFR_GLR; + } + else if (gfr_string == "OGR") { + gfr_type_ = GFR_OGR; + } + else { + gfr_type_ = GFR_INVALID; + } } catch (std::invalid_argument& e) { - //TODO: log here - gfr_type_ = GFR_INVALID; + //TODO: log here + gfr_type_ = GFR_INVALID; } //Artificial lift try { - std::string alq_string = header->getItem("ALQ")->getString(0); - if (alq_string == "GRAT") { - alq_type_ = ALQ_GRAT; - } - else if (alq_string == "IGLR") { - alq_type_ = ALQ_IGLR; - } - else if (alq_string == "TGLR") { - alq_type_ = ALQ_TGLR; - } - else if (alq_string == "PUMP") { - alq_type_ = ALQ_PUMP; - } - else if (alq_string == "COMP") { - alq_type_ = ALQ_COMP; - } - else if (alq_string == "BEAN") { - alq_type_ = ALQ_BEAN; - } - else if (alq_string == "UNDEF") { - alq_type_ = ALQ_UNDEF; - } - else { - alq_type_ = ALQ_INVALID; - } + std::string alq_string = header->getItem("ALQ")->getString(0); + if (alq_string == "GRAT") { + alq_type_ = ALQ_GRAT; + } + else if (alq_string == "IGLR") { + alq_type_ = ALQ_IGLR; + } + else if (alq_string == "TGLR") { + alq_type_ = ALQ_TGLR; + } + else if (alq_string == "PUMP") { + alq_type_ = ALQ_PUMP; + } + else if (alq_string == "COMP") { + alq_type_ = ALQ_COMP; + } + else if (alq_string == "BEAN") { + alq_type_ = ALQ_BEAN; + } + else if (alq_string == "UNDEF") { + alq_type_ = ALQ_UNDEF; + } + else { + alq_type_ = ALQ_INVALID; + } + } + catch (std::invalid_argument& e) { + //TODO: log here + alq_type_ = ALQ_INVALID; } - catch (std::invalid_argument& e) { - //TODO: log here - alq_type_ = ALQ_INVALID; - } //Get actual rate / flow values flo_data_ = (*iter++)->getItem("FLOW_VALUES")->getRawDoubleData(); @@ -304,8 +304,8 @@ namespace Opm { //Find interpolation ratio double dist = (*ceil_iter - *floor_iter); if (dist > 0) { - //Possible source for floating point error here if value and floor are large, - //but very close to each other + //Possible source for floating point error here if value and floor are large, + //but very close to each other retval.factor_ = (value-*floor_iter) / dist; } else { @@ -319,7 +319,7 @@ namespace Opm { const InterpData& wfr_i, const InterpData& gfr_i, const InterpData& alq_i) { //extents shape({{2, 2, 2, 2, 2}}); //array_type nn(shape); - double nn[2][2][2][2][2]; + double nn[2][2][2][2][2]; //Pick out nearest neighbors (nn) to our evaluation point //The following ladder of for loops will presumably be unrolled by a reasonable compiler. diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 01641e996..257715a67 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -48,72 +48,72 @@ extern const double reference[]; */ BOOST_AUTO_TEST_CASE(InterpolateZero) { - //All of our axes go from 0 to 1 - const std::vector& thp_axis{0.0, 1.0}; - const std::vector& wfr_axis{0.0, 0.5, 1.0}; - const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; - const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; - const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, - 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; + //All of our axes go from 0 to 1 + const std::vector& thp_axis{0.0, 1.0}; + const std::vector& wfr_axis{0.0, 0.5, 1.0}; + const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; + const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; + const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, + 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; - int nx = static_cast(thp_axis.size()); - int ny = static_cast(wfr_axis.size()); - int nz = static_cast(gfr_axis.size()); - int nu = static_cast(alq_axis.size()); - int nv = static_cast(flo_axis.size()); + int nx = static_cast(thp_axis.size()); + int ny = static_cast(wfr_axis.size()); + int nz = static_cast(gfr_axis.size()); + int nu = static_cast(alq_axis.size()); + int nv = static_cast(flo_axis.size()); - //Allocate data and fill with trivial data (zeros) - Opm::VFPProperties::extents size = {{ nx, ny, nz, nu, nv }}; - Opm::VFPProperties::array_type data(size); - for (int i=0; i(n-1); - for (int j=0; j(n-1); - for (int k=0; k(n-1); - for (int l=0; l(n-1); - for (int m=0; m(n-1); + //Check interpolation + double sum = 0.0; + int n=5; + for (int i=0; i(n-1); + for (int j=0; j(n-1); + for (int k=0; k(n-1); + for (int l=0; l(n-1); + for (int m=0; m(n-1); - //Note order of arguments! - sum += table.bhp(v, x, y, z, u); - } - } - } - } - } + //Note order of arguments! + sum += table.bhp(v, x, y, z, u); + } + } + } + } + } - BOOST_CHECK_EQUAL(sum, 0.0); + BOOST_CHECK_EQUAL(sum, 0.0); } @@ -123,73 +123,73 @@ BOOST_AUTO_TEST_CASE(InterpolateZero) */ BOOST_AUTO_TEST_CASE(InterpolateOne) { - //All of our axes go from 0 to 1 - const std::vector& thp_axis{0.0, 1.0}; - const std::vector& wfr_axis{0.0, 0.5, 1.0}; - const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; - const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; - const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, - 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; + //All of our axes go from 0 to 1 + const std::vector& thp_axis{0.0, 1.0}; + const std::vector& wfr_axis{0.0, 0.5, 1.0}; + const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; + const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; + const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, + 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; - int nx = static_cast(thp_axis.size()); - int ny = static_cast(wfr_axis.size()); - int nz = static_cast(gfr_axis.size()); - int nu = static_cast(alq_axis.size()); - int nv = static_cast(flo_axis.size()); + int nx = static_cast(thp_axis.size()); + int ny = static_cast(wfr_axis.size()); + int nz = static_cast(gfr_axis.size()); + int nu = static_cast(alq_axis.size()); + int nv = static_cast(flo_axis.size()); - //Allocate data and fill with trivial data (zeros) - Opm::VFPProperties::extents size = {{ nx, ny, nz, nu, nv }}; - Opm::VFPProperties::array_type data(size); - for (int i=0; i(n-1); - for (int j=0; j(n-1); - for (int k=0; k(n-1); - for (int l=0; l(n-1); - for (int m=0; m(n-1); + //Check interpolation + double sum = 0.0; + int n=5; + for (int i=0; i(n-1); + for (int j=0; j(n-1); + for (int k=0; k(n-1); + for (int l=0; l(n-1); + for (int m=0; m(n-1); - //Note order of arguments! - sum += table.bhp(v, x, y, z, u); - } - } - } - } - } + //Note order of arguments! + sum += table.bhp(v, x, y, z, u); + } + } + } + } + } - double reference = n*n*n*n*n; - BOOST_CHECK_EQUAL(sum, reference); + double reference = n*n*n*n*n; + BOOST_CHECK_EQUAL(sum, reference); } @@ -199,91 +199,91 @@ BOOST_AUTO_TEST_CASE(InterpolateOne) */ BOOST_AUTO_TEST_CASE(InterpolatePlane) { - //All of our axes go from 0 to 1 - const std::vector& thp_axis{0.0, 1.0}; - const std::vector& wfr_axis{0.0, 0.5, 1.0}; - const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; - const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; - const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, - 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; + //All of our axes go from 0 to 1 + const std::vector& thp_axis{0.0, 1.0}; + const std::vector& wfr_axis{0.0, 0.5, 1.0}; + const std::vector& gfr_axis{0.0, 0.25, 0.5, 0.75, 1}; + const std::vector& alq_axis{0.0, 0.125, 0.25, 0.375, 0.5, 0.625, 0.75, 0.875, 1}; + const std::vector& flo_axis{0.0, 0.0625, 0.125, 0.1875, 0.25, 0.3125, 0.375, 0.4375, + 0.5, 0.5625, 0.625, 0.6875, 0.75, 0.8125, 0.875, 0.9375, 1}; - int nx = static_cast(thp_axis.size()); - int ny = static_cast(wfr_axis.size()); - int nz = static_cast(gfr_axis.size()); - int nu = static_cast(alq_axis.size()); - int nv = static_cast(flo_axis.size()); + int nx = static_cast(thp_axis.size()); + int ny = static_cast(wfr_axis.size()); + int nz = static_cast(gfr_axis.size()); + int nu = static_cast(alq_axis.size()); + int nv = static_cast(flo_axis.size()); - //Allocate data and fill with trivial data (zeros) - Opm::VFPProperties::extents size = {{ nx, ny, nz, nu, nv }}; - Opm::VFPProperties::array_type data(size); - for (int i=0; i(nx-1); - for (int j=0; j(ny-1); - for (int k=0; k(nz-1); - for (int l=0; l(nu-1); - for (int m=0; m(nv-1); + //Allocate data and fill with trivial data (zeros) + Opm::VFPProperties::extents size = {{ nx, ny, nz, nu, nv }}; + Opm::VFPProperties::array_type data(size); + for (int i=0; i(nx-1); + for (int j=0; j(ny-1); + for (int k=0; k(nz-1); + for (int l=0; l(nu-1); + for (int m=0; m(nv-1); - data[i][j][k][l][m] = x + 2*y + 3*z + 4*u + 5*v; - } - } - } - } - } + data[i][j][k][l][m] = x + 2*y + 3*z + 4*u + 5*v; + } + } + } + } + } - //Create table - Opm::VFPProperties table(1, - 1000.0, - Opm::VFPProperties::FLO_OIL, - Opm::VFPProperties::WFR_WOR, - Opm::VFPProperties::GFR_GOR, - Opm::VFPProperties::ALQ_UNDEF, - flo_axis, - thp_axis, - wfr_axis, - gfr_axis, - alq_axis, - data); + //Create table + Opm::VFPProperties table(1, + 1000.0, + Opm::VFPProperties::FLO_OIL, + Opm::VFPProperties::WFR_WOR, + Opm::VFPProperties::GFR_GOR, + Opm::VFPProperties::ALQ_UNDEF, + flo_axis, + thp_axis, + wfr_axis, + gfr_axis, + alq_axis, + data); - //Check interpolation - double sum = 0.0; - double reference_sum = 0.0; - double sad = 0.0; // Sum absolute difference - double max_d = 0.0; // Maximum difference - int n=5; - for (int i=0; i<=n; ++i) { - const double x = i / static_cast(n); - for (int j=0; j<=n; ++j) { - const double y = j / static_cast(n); - for (int k=0; k<=n; ++k) { - const double z = k / static_cast(n); - for (int l=0; l<=n; ++l) { - const double u = l / static_cast(n); - for (int m=0; m<=n; ++m) { - const double v = m / static_cast(n); - double reference = x + 2*y + 3*z + 4*u + 5*v; - reference_sum += reference; + //Check interpolation + double sum = 0.0; + double reference_sum = 0.0; + double sad = 0.0; // Sum absolute difference + double max_d = 0.0; // Maximum difference + int n=5; + for (int i=0; i<=n; ++i) { + const double x = i / static_cast(n); + for (int j=0; j<=n; ++j) { + const double y = j / static_cast(n); + for (int k=0; k<=n; ++k) { + const double z = k / static_cast(n); + for (int l=0; l<=n; ++l) { + const double u = l / static_cast(n); + for (int m=0; m<=n; ++m) { + const double v = m / static_cast(n); + double reference = x + 2*y + 3*z + 4*u + 5*v; + reference_sum += reference; - //Note order of arguments! - double value = table.bhp(v, x, y, z, u); - sum += value; + //Note order of arguments! + double value = table.bhp(v, x, y, z, u); + sum += value; - double abs_diff = std::abs(value - reference); + double abs_diff = std::abs(value - reference); - sad += std::abs(abs_diff); - max_d = std::max(max_d, abs_diff); - } - } - } - } - } + sad += std::abs(abs_diff); + max_d = std::max(max_d, abs_diff); + } + } + } + } + } - BOOST_CHECK_CLOSE(sum, reference_sum, 0.0001); - BOOST_CHECK_SMALL(max_d, 0.0005); - BOOST_CHECK_SMALL(sad, 0.1); + BOOST_CHECK_CLOSE(sum, reference_sum, 0.0001); + BOOST_CHECK_SMALL(max_d, 0.0005); + BOOST_CHECK_SMALL(sad, 0.1); } @@ -297,70 +297,70 @@ BOOST_AUTO_TEST_CASE(Integration) Opm::DeckConstPtr deck; std::shared_ptr eclipseState; - Opm::ParserPtr parser(new Opm::Parser()); - boost::filesystem::path file("../opm-parser/testdata/integration_tests/VFPPROD/VFPPROD1"); + Opm::ParserPtr parser(new Opm::Parser()); + boost::filesystem::path file("../opm-parser/testdata/integration_tests/VFPPROD/VFPPROD1"); - deck = parser->parseFile(file.string()); - Opm::checkDeck(deck); + deck = parser->parseFile(file.string()); + Opm::checkDeck(deck); - BOOST_REQUIRE(deck->hasKeyword("VFPPROD")); + BOOST_REQUIRE(deck->hasKeyword("VFPPROD")); - std::vector tables; + std::vector tables; - int num_tables = deck->numKeywords("VFPPROD"); - for (int i=0; igetKeyword("VFPPROD" , i); - tables.push_back(Opm::VFPProperties(table)); - } + int num_tables = deck->numKeywords("VFPPROD"); + for (int i=0; igetKeyword("VFPPROD" , i); + tables.push_back(Opm::VFPProperties(table)); + } - //Do some rudimentary testing - //Get the BHP as a function of rate, thp, wfr, gfr, alq - double liq[] = {100, 2942.8571428571427, 5785.7142857142853, 8628.5714285714294, 11471.428571428571, 14314.285714285714, 17157.142857142859, 20000}; - double gor[] = {90, 1505.7142857142858, 2921.4285714285716, 4337.1428571428569, 5752.8571428571431, 7168.5714285714284, 8584.2857142857138, 10000}; - double wct[] = {0, 0.14285714285714285, 0.2857142857142857, 0.42857142857142855, 0.5714285714285714, 0.7142857142857143, 0.8571428571428571, 1}; - double thp[] = {16.010000000000002, 22.438571428571429, 28.867142857142859, 35.295714285714283, 41.724285714285713, 48.152857142857144, 54.581428571428575, 61.009999999999998}; - int n = sizeof(liq) / sizeof(liq[0]); + //Do some rudimentary testing + //Get the BHP as a function of rate, thp, wfr, gfr, alq + double liq[] = {100, 2942.8571428571427, 5785.7142857142853, 8628.5714285714294, 11471.428571428571, 14314.285714285714, 17157.142857142859, 20000}; + double gor[] = {90, 1505.7142857142858, 2921.4285714285716, 4337.1428571428569, 5752.8571428571431, 7168.5714285714284, 8584.2857142857138, 10000}; + double wct[] = {0, 0.14285714285714285, 0.2857142857142857, 0.42857142857142855, 0.5714285714285714, 0.7142857142857143, 0.8571428571428571, 1}; + double thp[] = {16.010000000000002, 22.438571428571429, 28.867142857142859, 35.295714285714283, 41.724285714285713, 48.152857142857144, 54.581428571428575, 61.009999999999998}; + int n = sizeof(liq) / sizeof(liq[0]); #ifdef verbose - std::cout << std::fixed << std::setprecision(17); - std::cout << "a=[..." << std::endl; + std::cout << std::fixed << std::setprecision(17); + std::cout << "a=[..." << std::endl; #endif - int i = 0; - double sad = 0.0; //Sum of absolute difference - double max_d = 0.0; //Maximum difference - for (int t=0; t