From b2335ced2485e4a25b9efcbc6117c62fa0a96f60 Mon Sep 17 00:00:00 2001 From: babrodtk Date: Wed, 19 Aug 2015 11:32:38 +0200 Subject: [PATCH] Minor fixes for PR --- opm/autodiff/BlackoilModelBase_impl.hpp | 1 - opm/autodiff/VFPHelpers.hpp | 17 ++--- opm/autodiff/VFPInjProperties.cpp | 4 +- opm/autodiff/VFPInjProperties.hpp | 6 +- opm/autodiff/VFPProdProperties.cpp | 56 ++++----------- opm/autodiff/VFPProdProperties.hpp | 6 +- opm/autodiff/VFPProperties.hpp | 6 +- opm/autodiff/WellDensitySegmented.cpp | 2 - tests/test_vfpproperties.cpp | 92 ++++++++++++++++++------- 9 files changed, 99 insertions(+), 91 deletions(-) diff --git a/opm/autodiff/BlackoilModelBase_impl.hpp b/opm/autodiff/BlackoilModelBase_impl.hpp index 4d2561bb9..0f6e7045e 100644 --- a/opm/autodiff/BlackoilModelBase_impl.hpp +++ b/opm/autodiff/BlackoilModelBase_impl.hpp @@ -1972,7 +1972,6 @@ namespace detail { vapour = wr[w*np + pu.phase_pos[ Gas ] ]; } - auto wc = wells().ctrls[w]; double alq = well_controls_iget_alq(wc, ctrl_index); int table_id = well_controls_iget_vfp(wc, ctrl_index); diff --git a/opm/autodiff/VFPHelpers.hpp b/opm/autodiff/VFPHelpers.hpp index 05cd61827..a279c1bc0 100644 --- a/opm/autodiff/VFPHelpers.hpp +++ b/opm/autodiff/VFPHelpers.hpp @@ -208,7 +208,7 @@ struct InterpData { inline InterpData findInterpData(const double& value, const std::vector& values) { InterpData retval; - const double abs_value = std::abs(value); + const double abs_value = value;//std::abs(value); //If we only have one value in our vector, return that if (values.size() == 1) { @@ -515,11 +515,12 @@ inline VFPEvaluation bhp(const VFPProdTable* table, double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); //First, find the values to interpolate between - auto flo_i = detail::findInterpData(flo, table->getFloAxis()); - auto thp_i = detail::findInterpData(thp, table->getTHPAxis()); - auto wfr_i = detail::findInterpData(wfr, table->getWFRAxis()); - auto gfr_i = detail::findInterpData(gfr, table->getGFRAxis()); - auto alq_i = detail::findInterpData(alq, table->getALQAxis()); + //Recall that flo is negative in Opm, so switch sign. + auto flo_i = detail::findInterpData(-flo, table->getFloAxis()); + auto thp_i = detail::findInterpData( thp, table->getTHPAxis()); + auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis()); + auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis()); + auto alq_i = detail::findInterpData( alq, table->getALQAxis()); detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); @@ -786,10 +787,10 @@ ADB gather_vars(const std::vector& well_tables, const ADB& values = map.find(key)->second; //Get indices to all elements that should use this ADB - const std::vector& elems = value; + const std::vector& current = value; //Add these elements to retval - retval = retval + superset(subset(values, elems), elems, values.size()); + retval = retval + superset(subset(values, current), current, values.size()); } return retval; diff --git a/opm/autodiff/VFPInjProperties.cpp b/opm/autodiff/VFPInjProperties.cpp index b8cc80d2a..6c4feb8a6 100644 --- a/opm/autodiff/VFPInjProperties.cpp +++ b/opm/autodiff/VFPInjProperties.cpp @@ -115,7 +115,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, ADB::V dflo = ADB::V::Zero(nw); //Get the table for each well - std::vector well_tables(nw, NULL); + std::vector well_tables(nw, nullptr); for (int i=0; i 0) { well_tables[i] = detail::getTable(m_tables, table_id[i]); @@ -128,7 +128,7 @@ VFPInjProperties::ADB VFPInjProperties::bhp(const std::vector& table_id, //Compute the BHP for each well independently for (int i=0; igetFloAxis()); auto thp_i = detail::findInterpData(thp.value()[i], table->getTHPAxis()); diff --git a/opm/autodiff/VFPInjProperties.hpp b/opm/autodiff/VFPInjProperties.hpp index c097aeba3..7a6e4c71e 100644 --- a/opm/autodiff/VFPInjProperties.hpp +++ b/opm/autodiff/VFPInjProperties.hpp @@ -46,14 +46,14 @@ public: * Takes *no* ownership of data. * @param inj_table A *single* VFPINJ table */ - VFPInjProperties(const VFPInjTable* inj_table); + explicit VFPInjProperties(const VFPInjTable* inj_table); /** * Constructor * Takes *no* ownership of data. * @param inj_tables A map of different VFPINJ tables. */ - VFPInjProperties(const std::map& inj_tables); + explicit VFPInjProperties(const std::map& inj_tables); /** * Linear interpolation of bhp as function of the input parameters. @@ -135,7 +135,7 @@ public: /** * Returns true if no vfp tables are in the current map */ - inline const bool empty() const { + bool empty() const { return m_tables.empty(); } diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index 995f6e9bf..dcaeef616 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -105,7 +105,7 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, ADB::V dflo = ADB::V::Zero(nw); //Get the table for each well - std::vector well_tables(nw, NULL); + std::vector well_tables(nw, nullptr); for (int i=0; i= 0) { well_tables[i] = detail::getTable(m_tables, table_id[i]); @@ -120,48 +120,17 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, //Compute the BHP for each well independently for (int i=0; igetFloAxis()); - auto thp_i = detail::findInterpData(thp.value()[i], table->getTHPAxis()); - auto wfr_i = detail::findInterpData(wfr.value()[i], table->getWFRAxis()); - auto gfr_i = detail::findInterpData(gfr.value()[i], table->getGFRAxis()); - auto alq_i = detail::findInterpData(alq.value()[i], table->getALQAxis()); + //Value of FLO is negative in OPM for producers, but positive in VFP table + auto flo_i = detail::findInterpData(-flo.value()[i], table->getFloAxis()); + auto thp_i = detail::findInterpData( thp.value()[i], table->getTHPAxis()); + auto wfr_i = detail::findInterpData( wfr.value()[i], table->getWFRAxis()); + auto gfr_i = detail::findInterpData( gfr.value()[i], table->getGFRAxis()); + auto alq_i = detail::findInterpData( alq.value()[i], table->getALQAxis()); detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); - /* - static const int N=40; - std::cout << "bhp=" << bhp_val.value << ";" << std::endl; - std::cout << "flo=" << flo.value()[i] << ";" << std::endl; - std::cout << "thp=" << thp.value()[i] << ";" << std::endl; - std::cout << "wfr=" << wfr.value()[i] << ";" << std::endl; - std::cout << "gfr=" << gfr.value()[i] << ";" << std::endl; - std::cout << "alq=" << alq.value()[i] << ";" << std::endl; - std::cout << "bhp_vfp=[" << std::endl; - for (int j=0; jgetFloAxis().front(); - const double end = table->getFloAxis().back(); - const double dist = end - start; - double flo_d = start + (j/static_cast(N-1)) * dist; - auto flo_i = detail::findInterpData(flo_d, table->getFloAxis()); - detail::VFPEvaluation bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); - std::cout << bhp_val.value << ","; - } - std::cout << "];" << std::endl; - - std::cout << "flo_vfp=[" << std::endl; - for (int j=0; jgetFloAxis().front(); - const double end = table->getFloAxis().back(); - const double dist = end - start; - double flo_d = start + (j/static_cast(N-1)) * dist; - std::cout << flo_d << ","; - } - std::cout << "];" << std::endl; - */ - - value[i] = bhp_val.value; dthp[i] = bhp_val.dthp; dwfr[i] = bhp_val.dwfr; @@ -247,11 +216,12 @@ double VFPProdProperties::thp(int table_id, * Find the function bhp_array(thp) by creating a 1D view of the data * by interpolating for every value of thp. This might be somewhat * expensive, but let us assome that nthp is small + * Recall that flo is negative in Opm, so switch the sign */ - auto flo_i = detail::findInterpData(flo, table->getFloAxis()); - auto wfr_i = detail::findInterpData(wfr, table->getWFRAxis()); - auto gfr_i = detail::findInterpData(gfr, table->getGFRAxis()); - auto alq_i = detail::findInterpData(alq, table->getALQAxis()); + auto flo_i = detail::findInterpData(-flo, table->getFloAxis()); + auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis()); + auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis()); + auto alq_i = detail::findInterpData( alq, table->getALQAxis()); std::vector bhp_array(nthp); for (int i=0; i& prod_tables); + explicit VFPProdProperties(const std::map& prod_tables); /** * Linear interpolation of bhp as function of the input parameters. @@ -146,7 +146,7 @@ public: /** * Returns true if no vfp tables are in the current map */ - inline const bool empty() const { + bool empty() const { return m_tables.empty(); } diff --git a/opm/autodiff/VFPProperties.hpp b/opm/autodiff/VFPProperties.hpp index 6387638e3..7d381f0ec 100644 --- a/opm/autodiff/VFPProperties.hpp +++ b/opm/autodiff/VFPProperties.hpp @@ -44,7 +44,7 @@ public: * @param inj_table A *single* VFPINJ table or NULL (no table) * @param prod_table A *single* VFPPROD table or NULL (no table) */ - VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table); + explicit VFPProperties(const VFPInjTable* inj_table, const VFPProdTable* prod_table); /** * Constructor @@ -58,14 +58,14 @@ public: /** * Returns the VFP properties for injection wells */ - inline const VFPInjProperties* getInj() const { + const VFPInjProperties* getInj() const { return m_inj.get(); } /** * Returns the VFP properties for production wells */ - inline const VFPProdProperties* getProd() const { + const VFPProdProperties* getProd() const { return m_prod.get(); } diff --git a/opm/autodiff/WellDensitySegmented.cpp b/opm/autodiff/WellDensitySegmented.cpp index 576fb8583..ae43c9532 100644 --- a/opm/autodiff/WellDensitySegmented.cpp +++ b/opm/autodiff/WellDensitySegmented.cpp @@ -147,7 +147,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, const std::vector& z_perf, const std::vector& dens_perf, const double gravity) { - const int np = wells.number_of_phases; const int nw = wells.number_of_wells; const int nperf = wells.well_connpos[nw]; @@ -177,7 +176,6 @@ Opm::WellDensitySegmented::computeConnectionPressureDelta(const Wells& wells, const double z_above = perf == wells.well_connpos[w] ? wells.depth_ref[w] : z_perf[perf - 1]; const double dz = z_perf[perf] - z_above; dp_perf[perf] = dz * dens_perf[perf] * gravity; - //dens[wells.well_connpos[w]] } } diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index bde853043..06bbbc5eb 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -31,9 +31,12 @@ #include #include #include +#include +#include #include #include +#include #include #include @@ -58,6 +61,43 @@ const double sad_tol = 1.0e-8; + +BOOST_AUTO_TEST_SUITE( HelperTests ) + +BOOST_AUTO_TEST_CASE(findInterpData) +{ + std::vector values = {1, 5, 7, 9, 11, 15}; + double interpolate = 6.87; + double extrapolate_left = -1.89; + double extrapolate_right = 32.1; + + Opm::detail::InterpData eval1 = Opm::detail::findInterpData(interpolate, values); + Opm::detail::InterpData eval2 = Opm::detail::findInterpData(extrapolate_left, values); + Opm::detail::InterpData eval3 = Opm::detail::findInterpData(extrapolate_right, values); + + BOOST_CHECK_EQUAL(eval1.ind_[0], 1); + BOOST_CHECK_EQUAL(eval1.ind_[1], 2); + BOOST_CHECK_EQUAL(eval1.factor_, (interpolate-values[1]) / (values[2] - values[1])); + + BOOST_CHECK_EQUAL(eval2.ind_[0], 0); + BOOST_CHECK_EQUAL(eval2.ind_[1], 1); + BOOST_CHECK_EQUAL(eval2.factor_, (extrapolate_left-values[0]) / (values[1] - values[0])); + + BOOST_CHECK_EQUAL(eval3.ind_[0], 4); + BOOST_CHECK_EQUAL(eval3.ind_[1], 5); + BOOST_CHECK_EQUAL(eval3.factor_, (extrapolate_right-values[4]) / (values[5] - values[4])); +} + +BOOST_AUTO_TEST_SUITE_END() // HelperTests + + + + + + + + + struct ConversionFixture { typedef Opm::VFPProdProperties::ADB ADB; @@ -391,11 +431,11 @@ BOOST_AUTO_TEST_CASE(GetTable) add_well(INJECTOR, 100, 1, NULL, cells, NULL, NULL, wells.get()); //Create interpolation points - double aqua_d = 0.15; - double liquid_d = 0.25; - double vapour_d = 0.35; - double thp_d = 0.45; - double alq_d = 0.55; + double aqua_d = -0.15; + double liquid_d = -0.25; + double vapour_d = -0.35; + double thp_d = 0.45; + double alq_d = 0.55; ADB aqua_adb = createConstantScalarADB(aqua_d); ADB liquid_adb = createConstantScalarADB(liquid_d); @@ -408,7 +448,7 @@ BOOST_AUTO_TEST_CASE(GetTable) ADB qs_adb = ADB::constant(qs_adb_v); //Check that our reference has not changed - Opm::detail::VFPEvaluation ref= Opm::detail::bhp(&table, aqua_d, liquid_d, vapour_d, thp_d, alq_d); + Opm::detail::VFPEvaluation ref = Opm::detail::bhp(&table, aqua_d, liquid_d, vapour_d, thp_d, alq_d); BOOST_CHECK_CLOSE(ref.value, 1.0923565702101556, max_d_tol); BOOST_CHECK_CLOSE(ref.dthp, 0.13174065498177251, max_d_tol); BOOST_CHECK_CLOSE(ref.dwfr, -1.2298177745501071, max_d_tol); @@ -523,13 +563,13 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane) for (int i=0; i<=n; ++i) { const double thp = i / static_cast(n); for (int j=1; j<=n; ++j) { - const double aqua = j / static_cast(n); + const double aqua = -j / static_cast(n); for (int k=1; k<=n; ++k) { - const double vapour = k / static_cast(n); + const double vapour = -k / static_cast(n); for (int l=0; l<=n; ++l) { const double alq = l / static_cast(n); for (int m=1; m<=n; ++m) { - const double liquid = m / static_cast(n); + const double liquid = -m / static_cast(n); //Find values that should be in table double flo = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType()); @@ -537,7 +577,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane) double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType()); //Calculate reference - double reference = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo; + double reference = thp + 2*wfr + 3*gfr+ 4*alq - 5*flo; //Calculate actual //Note order of arguments: id, aqua, liquid, vapour, thp, alq @@ -545,7 +585,7 @@ BOOST_AUTO_TEST_CASE(InterpolatePlane) double abs_diff = std::abs(actual - reference); - double max_d = std::max(max_d, abs_diff); + max_d = std::max(max_d, abs_diff); sad = sad + abs_diff; } } @@ -578,20 +618,20 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane) for (int i=0; i<=n+o; ++i) { const double x = i / static_cast(n); for (int j=1; j<=n+o; ++j) { - const double aqua = j / static_cast(n); + const double aqua = -j / static_cast(n); for (int k=1; k<=n+o; ++k) { - const double vapour = k / static_cast(n); + const double vapour = -k / static_cast(n); for (int l=0; l<=n+o; ++l) { const double u = l / static_cast(n); for (int m=1; m<=n+o; ++m) { - const double liquid = m / static_cast(n); + const double liquid = -m / static_cast(n); //Find values that should be in table double v = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType()); double y = Opm::detail::getWFR(aqua, liquid, vapour, table.getWFRType()); double z = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType()); - double reference = x + 2*y + 3*z+ 4*u + 5*v; + double reference = x + 2*y + 3*z+ 4*u - 5*v; reference_sum += reference; //Note order of arguments! id, aqua, liquid, vapour, thp , alq @@ -635,13 +675,13 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB) for (int i=0; i<=n+o; ++i) { const double x = i / static_cast(n); for (int j=1; j<=n+o; ++j) { - const double aqua = j / static_cast(n); + const double aqua = -j / static_cast(n); for (int k=1; k<=n+o; ++k) { - const double vapour = k / static_cast(n); + const double vapour = -k / static_cast(n); for (int l=0; l<=n+o; ++l) { const double u = l / static_cast(n); for (int m=1; m<=n+o; ++m) { - const double liquid = m / static_cast(n); + const double liquid = -m / static_cast(n); //Temporary variables used to represent independent wells const int num_wells = 5; @@ -678,7 +718,7 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB) double y = Opm::detail::getWFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getWFRType()); double z = Opm::detail::getGFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getGFRType()); - reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) + 5*v; + reference = x*(w+1) + 2*y + 3*z + 4*u*(w+1) - 5*v; value = bhp_val[w]; sum += value; @@ -733,7 +773,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs) ADB::V qs_v(nphases*nwells); for (int j=0; j(nwells*nphases-1.0); + qs_v[j*nwells+i] = -(j*nwells+i) / static_cast(nwells*nphases-1.0); } } ADB qs = ADB::constant(qs_v); @@ -778,7 +818,7 @@ BOOST_AUTO_TEST_CASE(InterpolateADBAndQs) double flo = oil[i]; double wor = water[i]/oil[i]; double gor = gas[i]/oil[i]; - reference[i] = thp_v[i] + 2*wor + 3*gor + 4*alq_v[i] + 5*flo; + reference[i] = thp_v[i] + 2*wor + 3*gor + 4*alq_v[i] - 5*flo; } //Check that interpolation matches @@ -817,13 +857,13 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives) for (int i=0; i<=n; ++i) { const double thp = i / static_cast(n); for (int j=1; j<=n; ++j) { - const double aqua = j / static_cast(n); + const double aqua = -j / static_cast(n); for (int k=1; k<=n; ++k) { - const double vapour = k / static_cast(n); + const double vapour = -k / static_cast(n); for (int l=0; l<=n; ++l) { const double alq = l / static_cast(n); for (int m=1; m<=n; ++m) { - const double liquid = m / static_cast(n); + const double liquid = -m / static_cast(n); //Find values that should be in table double flo = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType()); @@ -832,7 +872,7 @@ BOOST_AUTO_TEST_CASE(PartialDerivatives) //Calculate reference VFPEvaluation reference; - reference.value = thp + 2*wfr + 3*gfr+ 4*alq + 5*flo; + reference.value = thp + 2*wfr + 3*gfr+ 4*alq - 5*flo; reference.dthp = 1; reference.dwfr = 2; reference.dgfr = 3; @@ -1040,7 +1080,7 @@ BOOST_AUTO_TEST_CASE(ParseInterpolateRealisticVFPPROD) //for (unsigned int a=0; a convert to SM3/second - double f_i = liq[f]*1.1574074074074073e-05; + double f_i = -liq[f]*1.1574074074074073e-05; //THP given as BARSA => convert to Pascal double t_i = thp[t]*100000.0;