diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index eb38e90cd..6feed3db6 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -55,6 +55,7 @@ list (APPEND TEST_SOURCE_FILES tests/test_scalar_mult.cpp tests/test_transmissibilitymultipliers.cpp tests/test_welldensitysegmented.cpp + tests/test_vfpproperties.cpp ) list (APPEND TEST_DATA_FILES diff --git a/opm/autodiff/VFPProperties.hpp b/opm/autodiff/VFPProperties.hpp index 7805f86f4..11f484852 100644 --- a/opm/autodiff/VFPProperties.hpp +++ b/opm/autodiff/VFPProperties.hpp @@ -26,121 +26,200 @@ namespace Opm { class VFPProperties { - public: - VFPProperties(DeckKeywordConstPtr table) { + public: + typedef boost::multi_array array_type; + typedef boost::array extents; + + ///Rate type + enum FLO_TYPE { + FLO_OIL, //< Oil rate + FLO_LIQ, //< Liquid rate + FLO_GAS, //< Gas rate + //FLO_WG + //FLO_TM + FLO_INVALID + }; + + ///Water fraction variable + enum WFR_TYPE { + WFR_WOR, //< Water-oil ratio + WFR_WCT, //< Water cut + WFR_WGR, //< Water-gas ratio + WFR_INVALID + }; + + ///Gas fraction variable + enum GFR_TYPE { + GFR_GOR, //< Gas-oil ratio + GFR_GLR, //< Gas-liquid ratio + GFR_OGR, //< Oil-gas ratio + GFR_INVALID + }; + + ///Artificial lift quantity + enum ALQ_TYPE { + ALQ_GRAT, //< Lift as injection rate + ALQ_IGLR, //< Injection gas-liquid ratio + ALQ_TGLR, //< Total gas-liquid ratio + ALQ_PUMP, //< Pump rating + ALQ_COMP, //< Compressor power + ALQ_BEAN, //< Choke diameter + ALQ_UNDEF, //< Undefined + ALQ_INVALID + }; + + 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) { + + } + + VFPProperties(DeckKeywordConstPtr table) { auto iter = table->begin(); auto header = (*iter++); - table_num_ = header->getItem("TABLE")->getInt(0); - datum_depth_ = header->getItem("DATUM_DEPTH")->getRawDouble(0); + table_num_ = header->getItem("TABLE")->getInt(0); + datum_depth_ = header->getItem("DATUM_DEPTH")->getRawDouble(0); - //Rate type - 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; - } - - //Water fraction - 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; + //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; + } + } + catch (std::invalid_argument& e) { + //TODO: log here + flo_type_ = FLO_INVALID; } - else { - wfr_type_ = WFR_INVALID; - } - //Gas fraction - 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; - } - - //Artificial lift - /* - * ALQ not implemented properly in parser? - - 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; + //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; + } + } + catch (std::invalid_argument& e) { + //TODO: log here + wfr_type_ = WFR_INVALID; } - else if (alq_string == "PUMP") { - alq_type_ = ALQ_PUMP; + + //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; + } + } + catch (std::invalid_argument& e) { + //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; + } + } + catch (std::invalid_argument& e) { + //TODO: log here + alq_type_ = ALQ_INVALID; } - 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; - } - */ - //Get actual rate / flow values - const std::vector& flo = (*iter++)->getItem("FLOW_VALUES")->getRawDoubleData(); - flo_data_.resize(flo.size()); - std::copy(flo.begin(), flo.end(), flo_data_.begin()); + //Get actual rate / flow values + flo_data_ = (*iter++)->getItem("FLOW_VALUES")->getRawDoubleData(); - //Get actual tubing head pressure values - const std::vector& thp = (*iter++)->getItem("THP_VALUES")->getRawDoubleData(); - thp_data_.resize(thp.size()); - std::copy(thp.begin(), thp.end(), thp_data_.begin()); + //Get actual tubing head pressure values + thp_data_ = (*iter++)->getItem("THP_VALUES")->getRawDoubleData(); - //Get actual water fraction values - const std::vector& wfr = (*iter++)->getItem("WFR_VALUES")->getRawDoubleData(); - wfr_data_.resize(wfr.size()); - std::copy(wfr.begin(), wfr.end(), wfr_data_.begin()); + //Get actual water fraction values + wfr_data_ = (*iter++)->getItem("WFR_VALUES")->getRawDoubleData(); - //Get actual gas fraction values - const std::vector& gfr = (*iter++)->getItem("GFR_VALUES")->getRawDoubleData(); - gfr_data_.resize(gfr.size()); - std::copy(gfr.begin(), gfr.end(), gfr_data_.begin()); + //Get actual gas fraction values + gfr_data_ = (*iter++)->getItem("GFR_VALUES")->getRawDoubleData(); - //Get actual gas fraction values - const std::vector& alq = (*iter++)->getItem("ALQ_VALUES")->getRawDoubleData(); - alq_data_.resize(alq.size()); - std::copy(alq.begin(), alq.end(), alq_data_.begin()); + //Get actual gas fraction values + alq_data_ = (*iter++)->getItem("ALQ_VALUES")->getRawDoubleData(); //Finally, read the actual table itself. - unsigned int nt = thp_data_.size(); - unsigned int nw = wfr_data_.size(); - unsigned int ng = gfr_data_.size(); - unsigned int na = alq_data_.size(); - unsigned int nf = flo_data_.size(); + size_t nt = thp_data_.size(); + size_t nw = wfr_data_.size(); + size_t ng = gfr_data_.size(); + size_t na = alq_data_.size(); + size_t nf = flo_data_.size(); extents shape; shape[0] = nt; shape[1] = nw; @@ -149,218 +228,181 @@ namespace Opm { shape[4] = nf; data_.resize(shape); - for (; iter!=table->end(); ++iter) { - //Get indices (subtract 1 to get 0-based index) - unsigned int t = (*iter)->getItem("THP_INDEX")->getInt(0) - 1; - unsigned int w = (*iter)->getItem("WFR_INDEX")->getInt(0) - 1; - unsigned int g = (*iter)->getItem("GFR_INDEX")->getInt(0) - 1; - unsigned int a = (*iter)->getItem("ALQ_INDEX")->getInt(0) - 1; + for (; iter!=table->end(); ++iter) { + //Get indices (subtract 1 to get 0-based index) + int t = (*iter)->getItem("THP_INDEX")->getInt(0) - 1; + int w = (*iter)->getItem("WFR_INDEX")->getInt(0) - 1; + int g = (*iter)->getItem("GFR_INDEX")->getInt(0) - 1; + int a = (*iter)->getItem("ALQ_INDEX")->getInt(0) - 1; - //Rest of values (bottom hole pressure or tubing head temperature) have index of flo value - const std::vector& bhp_tht = (*iter)->getItem("VALUES")->getRawDoubleData(); - std::copy(bhp_tht.begin(), bhp_tht.end(), &data_[t][w][g][a][0]); + //Rest of values (bottom hole pressure or tubing head temperature) have index of flo value + const std::vector& bhp_tht = (*iter)->getItem("VALUES")->getRawDoubleData(); + std::copy(bhp_tht.begin(), bhp_tht.end(), &data_[t][w][g][a][0]); - //Check for large values - for (unsigned int i = 0; i 1.0e10) { - //TODO: Replace with proper log message - std::cerr << "Too large value encountered in VFPPROD in [" - << t << "," << w << "," << g << "," << a << "]=" - << bhp_tht[i] << std::endl; - } - } + //Check for large values + for (size_t i = 0; i 1.0e10) { + //TODO: Replace with proper log message + std::cerr << "Too large value encountered in VFPPROD in [" + << t << "," << w << "," << g << "," << a << "]=" + << bhp_tht[i] << std::endl; + } + } } - } + } - struct InterpData { - InterpData() : ind_({0}), factor_(0.0) {} - unsigned int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value] - double factor_; //Interpolation factor - }; + /** + * Linear interpolation of bhp as a function of the input parameters + * @param flo Production rate of oil, gas or liquid + * @param thp Tubing head pressure + * @param wfr Water-oil ratio, water cut, or water-gas ratio + * @param gfr Gas-oil ratio, gas-liquid ratio, or oil-gas ratio + * @param alq Artificial lift or other parameter + * + * @return The bottom hole pressure, interpolated linearly using + * the above parameters from the values in the input table. + */ + double bhp(double flo, double thp, float wfr, float gfr, float alq) { + //First, find the floor value of the inputs + auto flo_i = find_interp_data(flo, flo_data_); + auto thp_i = find_interp_data(thp, thp_data_); + auto wfr_i = find_interp_data(wfr, wfr_data_); + auto gfr_i = find_interp_data(gfr, gfr_data_); + auto alq_i = find_interp_data(alq, alq_data_); - InterpData find_interp_data(double value, const std::vector& values) { - InterpData retval; + return interpolate(flo_i, thp_i, wfr_i, gfr_i, alq_i); + } - //First element greater than or equal to value - //Don't access out-of-range, therefore values.end()-1 - auto ceil_iter = std::lower_bound(values.begin(), values.end()-1, value); + private: + struct InterpData { + InterpData() : factor_(0.0) {} + int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value] + double factor_; //Interpolation factor + }; - //Find last element smaller than or equal to range - auto floor_iter = ceil_iter; - if (*floor_iter == value) { - // floor_iter == ceil_iter == value - } - else if (floor_iter > values.begin()) { - // floor_iter <= value <= ceil_iter - --floor_iter; - } + InterpData find_interp_data(double value, const std::vector& values) { + InterpData retval; - //Now set these in the retval struct - retval.ind_[0] = floor_iter - values.begin(); - retval.ind_[1] = ceil_iter - values.begin(); + //First element greater than or equal to value + //Don't access out-of-range, therefore values.end()-1 + auto ceil_iter = std::lower_bound(values.begin(), values.end()-1, value); - //Find interpolation ratio - double dist = (*ceil_iter - *floor_iter); - if (dist > 0) { - retval.factor_ = (value-*floor_iter) / dist; - } - else { - retval.factor_ = 1.0; - } + //Find last element smaller than or equal to range + auto floor_iter = ceil_iter; + if (*floor_iter == value) { + // floor_iter == ceil_iter == value + } + else if (floor_iter > values.begin()) { + // floor_iter <= value <= ceil_iter + --floor_iter; + } - return retval; - } + //Now set these in the retval struct + retval.ind_[0] = floor_iter - values.begin(); + retval.ind_[1] = ceil_iter - values.begin(); - double interpolate(const InterpData& flo_i, const InterpData& thp_i, - const InterpData& wfr_i, const InterpData& gfr_i, const InterpData& alq_i) { - extents shape({{2, 2, 2, 2, 2}}); - array_type nn(shape); + //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 + retval.factor_ = (value-*floor_iter) / dist; + } + else { + retval.factor_ = 1.0; + } - //Pick out nearest neighbors (nn) to our evaluation point - //The following ladder of for loops will presumably be unrolled by a reasonable compiler. - //This is not really required, but performance-wise it may pay off, since the 32-elements - //we copy to (nn) will fit better in cache than the full original table for the - //interpolation below. - for (unsigned int t=0; t<=1; ++t) { - for (unsigned int w=0; w<=1; ++w) { - for (unsigned int g=0; g<=1; ++g) { - for (unsigned int a=0; a<=1; ++a) { - for (unsigned int f=0; f<=1; ++f) { - //Shorthands for indexing - unsigned int ti = thp_i.ind_[t]; - unsigned int wi = wfr_i.ind_[w]; - unsigned int gi = gfr_i.ind_[g]; - unsigned int ai = alq_i.ind_[a]; - unsigned int fi = flo_i.ind_[f]; + return retval; + } - //Copy element - nn[t][w][g][a][f] = data_[ti][wi][gi][ai][fi]; - } - } - } - } - } + double interpolate(const InterpData& flo_i, const InterpData& thp_i, + 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]; - //Remove dimensions iteratively - // Example: going from 3D to 2D to 1D, we start by interpolating along - // the z axis first, leaving a 2D problem. Then interpolating along the y - // axis, leaving a 1D, problem, etc. - double tf = flo_i.factor_; - for (unsigned int t=0; t<=1; ++t) { - for (unsigned int w=0; w<=1; ++w) { - for (unsigned int g=0; g<=1; ++g) { - for (unsigned int a=0; a<=1; ++a) { - nn[t][w][g][a][0] = (1.0-tf)*nn[t][w][g][a][0] + tf*nn[t][w][g][a][1]; - } - } - } - } + //Pick out nearest neighbors (nn) to our evaluation point + //The following ladder of for loops will presumably be unrolled by a reasonable compiler. + //This is not really required, but performance-wise it may pay off, since the 32-elements + //we copy to (nn) will fit better in cache than the full original table for the + //interpolation below. + for (int t=0; t<=1; ++t) { + for (int w=0; w<=1; ++w) { + for (int g=0; g<=1; ++g) { + for (int a=0; a<=1; ++a) { + for (int f=0; f<=1; ++f) { + //Shorthands for indexing + int ti = thp_i.ind_[t]; + int wi = wfr_i.ind_[w]; + int gi = gfr_i.ind_[g]; + int ai = alq_i.ind_[a]; + int fi = flo_i.ind_[f]; - tf = alq_i.factor_; - for (unsigned int t=0; t<=1; ++t) { - for (unsigned int w=0; w<=1; ++w) { - for (unsigned int g=0; g<=1; ++g) { - nn[t][w][g][0][0] = (1.0-tf)*nn[t][w][g][0][0] + tf*nn[t][w][g][1][0]; - } - } - } + //Copy element + nn[t][w][g][a][f] = data_[ti][wi][gi][ai][fi]; + } + } + } + } + } - tf = gfr_i.factor_; - for (unsigned int t=0; t<=1; ++t) { - for (unsigned int w=0; w<=1; ++w) { - nn[t][w][0][0][0] = (1.0-tf)*nn[t][w][0][0][0] + tf*nn[t][w][1][0][0]; - } - } + //Remove dimensions iteratively + // Example: going from 3D to 2D to 1D, we start by interpolating along + // the z axis first, leaving a 2D problem. Then interpolating along the y + // axis, leaving a 1D, problem, etc. + double tf = flo_i.factor_; + for (int t=0; t<=1; ++t) { + for (int w=0; w<=1; ++w) { + for (int g=0; g<=1; ++g) { + for (int a=0; a<=1; ++a) { + nn[t][w][g][a][0] = (1.0-tf)*nn[t][w][g][a][0] + tf*nn[t][w][g][a][1]; + } + } + } + } - tf = wfr_i.factor_; - for (unsigned int t=0; t<=1; ++t) { - nn[t][0][0][0][0] = (1.0-tf)*nn[t][0][0][0][0] + tf*nn[t][1][0][0][0]; - } + tf = alq_i.factor_; + for (int t=0; t<=1; ++t) { + for (int w=0; w<=1; ++w) { + for (int g=0; g<=1; ++g) { + nn[t][w][g][0][0] = (1.0-tf)*nn[t][w][g][0][0] + tf*nn[t][w][g][1][0]; + } + } + } - tf = thp_i.factor_; - return (1.0-tf)*nn[0][0][0][0][0] + tf*nn[1][0][0][0][0]; - } + tf = gfr_i.factor_; + for (int t=0; t<=1; ++t) { + for (int w=0; w<=1; ++w) { + nn[t][w][0][0][0] = (1.0-tf)*nn[t][w][0][0][0] + tf*nn[t][w][1][0][0]; + } + } - /** - * Linear interpolation of bhp as a function of the input parameters - * @param flo Production rate of oil, gas or liquid - * @param thp Tubing head pressure - * @param wfr Water-oil ratio, water cut, or water-gas ratio - * @param gfr Gas-oil ratio, gas-liquid ratio, or oil-gas ratio - * @param alq Artificial lift or other parameter - * - * @return The bottom hole pressure, interpolated linearly using - * the above parameters from the values in the input table. - */ - double bhp(double flo, double thp, float wfr, float gfr, float alq) { - //First, find the floor value of the inputs - auto flo_i = find_interp_data(flo, flo_data_); - auto thp_i = find_interp_data(thp, thp_data_); - auto wfr_i = find_interp_data(wfr, wfr_data_); - auto gfr_i = find_interp_data(gfr, gfr_data_); - auto alq_i = find_interp_data(alq, alq_data_); + tf = wfr_i.factor_; + for (int t=0; t<=1; ++t) { + nn[t][0][0][0][0] = (1.0-tf)*nn[t][0][0][0][0] + tf*nn[t][1][0][0][0]; + } - return interpolate(flo_i, thp_i, wfr_i, gfr_i, alq_i); - } + tf = thp_i.factor_; + return (1.0-tf)*nn[0][0][0][0][0] + tf*nn[1][0][0][0][0]; + } - ///Rate type - enum FLO_TYPE { - FLO_OIL, - FLO_LIQ, - FLO_GAS, - //FLO_WG - //FLO_TM - FLO_INVALID - }; + //"Header" variables + int table_num_; + double datum_depth_; + FLO_TYPE flo_type_; + WFR_TYPE wfr_type_; + GFR_TYPE gfr_type_; + ALQ_TYPE alq_type_; - ///Water fraction variable - enum WFR_TYPE { - WFR_WOR, //< Water-oil ratio - WFR_WCT, //< Water cut - WFR_WGR, //< Water-gas ratio - WFR_INVALID - }; + //The actual table axes + std::vector flo_data_; + std::vector thp_data_; + std::vector wfr_data_; + std::vector gfr_data_; + std::vector alq_data_; - ///Gas fraction variable - enum GFR_TYPE { - GFR_GOR, //< Gas-oil ratio - GFR_GLR, //< Gas-liquid ratio - GFR_OGR, //< Oil-gas ratio - GFR_INVALID - }; - - ///Artificial lift quantity - enum ALQ_TYPE { - ALQ_GRAT, //< Lift as injection rate - ALQ_IGLR, //< Injection gas-liquid ratio - ALQ_TGLR, //< Total gas-liquid ratio - ALQ_PUMP, //< Pump rating - ALQ_COMP, //< Compressor power - ALQ_BEAN, //< Choke diameter - ALQ_UNDEF, //< Undefined - ALQ_INVALID - }; - - private: - //"Header" variables - int table_num_; - double datum_depth_; - FLO_TYPE flo_type_; - WFR_TYPE wfr_type_; - GFR_TYPE gfr_type_; - ALQ_TYPE alq_type_; - - //The actual table axes - std::vector flo_data_; - std::vector thp_data_; - std::vector wfr_data_; - std::vector gfr_data_; - std::vector alq_data_; - - //The data itself - typedef boost::multi_array array_type; - typedef boost::array extents; - array_type data_; + //The data itself + array_type data_; }; } diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp new file mode 100644 index 000000000..01641e996 --- /dev/null +++ b/tests/test_vfpproperties.cpp @@ -0,0 +1,372 @@ +/* + Copyright 2015 SINTEF ICT, Applied Mathematics. + + This file is part of the Open Porous Media project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + + + +#include + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE AutoDiffBlockTest + +#include +#include + +#include +#include + +#include +#include +#include +#include + + +extern const double reference[]; + + +/** + * Test that we can generate some dummy zero-data, and + * interpolate something the analytic solution + */ +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}; + + 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); + + //Note order of arguments! + sum += table.bhp(v, x, y, z, u); + } + } + } + } + } + + BOOST_CHECK_EQUAL(sum, 0.0); +} + + +/** + * Test that we can generate some dummy one-data, and + * interpolate something the analytic solution + */ +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}; + + 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); + + //Note order of arguments! + sum += table.bhp(v, x, y, z, u); + } + } + } + } + } + + double reference = n*n*n*n*n; + BOOST_CHECK_EQUAL(sum, reference); +} + + +/** + * Test that we can generate some dummy one-data, and + * interpolate something the analytic solution + */ +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}; + + 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); + + 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); + + //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; + + double abs_diff = std::abs(value - reference); + + 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); +} + + + + +/** + * Tests that we can actually parse some input data, and interpolate within that space + */ +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"); + + deck = parser->parseFile(file.string()); + Opm::checkDeck(deck); + + BOOST_REQUIRE(deck->hasKeyword("VFPPROD")); + + std::vector tables; + + 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]); + +#ifdef verbose + 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