diff --git a/opm/autodiff/VFPHelpers.hpp b/opm/autodiff/VFPHelpers.hpp new file mode 100644 index 000000000..1721b8146 --- /dev/null +++ b/opm/autodiff/VFPHelpers.hpp @@ -0,0 +1,358 @@ +/* + 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 . +*/ + + +#ifndef OPM_AUTODIFF_VFPHELPERS_HPP_ +#define OPM_AUTODIFF_VFPHELPERS_HPP_ + + +#include + + +/** + * This file contains a set of helper functions used by VFPProd / VFPInj. + */ +namespace Opm { +namespace detail { + + +typedef VFPProdProperties::ADB ADB; + + + + + + + + + + +/** + * Returns zero if input value is NaN + */ +inline double zeroIfNan(const double& value) { + return (std::isnan(value)) ? 0.0 : value; +} + + + + + +/** + * Returns zero for every entry in the ADB which is NaN + */ +inline ADB zeroIfNan(const ADB& values) { + Selector not_nan_selector(values.value(), Selector::NotNaN); + + const ADB::V z = ADB::V::Zero(values.size()); + const ADB zero = ADB::constant(z, values.blockPattern()); + + ADB retval = not_nan_selector.select(values, zero); + return retval; +} + + + + + +/** + * Computes the flo parameter according to the flo_type_ + * @return Production rate of oil, gas or liquid. + */ +template +static T getFlo(const T& aqua, const T& liquid, const T& vapour, + const VFPProdTable::FLO_TYPE& type) { + switch (type) { + case VFPProdTable::FLO_OIL: + //Oil = liquid phase + return liquid; + case VFPProdTable::FLO_LIQ: + //Liquid = aqua + liquid phases + return aqua + liquid; + case VFPProdTable::FLO_GAS: + //Gas = vapor phase + return vapour; + case VFPProdTable::FLO_INVALID: //Intentional fall-through + default: + OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'"); + } +} + + + + + + + +/** + * Computes the wfr parameter according to the wfr_type_ + * @return Production rate of oil, gas or liquid. + */ +template +static T getWFR(const T& aqua, const T& liquid, const T& vapour, + const VFPProdTable::WFR_TYPE& type) { + switch(type) { + case VFPProdTable::WFR_WOR: { + //Water-oil ratio = water / oil + T wor = aqua / liquid; + return zeroIfNan(wor); + } + case VFPProdTable::WFR_WCT: + //Water cut = water / (water + oil) + return zeroIfNan(aqua / (aqua + liquid)); + case VFPProdTable::WFR_WGR: + //Water-gas ratio = water / gas + return zeroIfNan(aqua / vapour); + case VFPProdTable::WFR_INVALID: //Intentional fall-through + default: + OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'"); + } +} + + + + + + +/** + * Computes the gfr parameter according to the gfr_type_ + * @return Production rate of oil, gas or liquid. + */ +template +static T getGFR(const T& aqua, const T& liquid, const T& vapour, + const VFPProdTable::GFR_TYPE& type) { + switch(type) { + case VFPProdTable::GFR_GOR: + // Gas-oil ratio = gas / oil + return zeroIfNan(vapour / liquid); + case VFPProdTable::GFR_GLR: + // Gas-liquid ratio = gas / (oil + water) + return zeroIfNan(vapour / (liquid + aqua)); + case VFPProdTable::GFR_OGR: + // Oil-gas ratio = oil / gas + return zeroIfNan(liquid / vapour); + case VFPProdTable::GFR_INVALID: //Intentional fall-through + default: + OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'"); + } +} + + + + + + +/** + * Helper struct for linear interpolation + */ +struct InterpData { + InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {} + int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value] + double inv_dist_; // 1 / distance between the two end points of the segment. Used to calculate derivatives and uses 1.0 / 0.0 = 0.0 as a convention + double factor_; // Interpolation factor +}; + + + + + + +/** + * Helper function to find indices etc. for linear interpolation + */ +inline InterpData findInterpData(const double& value, const std::vector& values) { + InterpData retval; + + //If we only have one value in our vector, return that + if (values.size() == 1) { + retval.ind_[0] = 0; + retval.ind_[1] = 0; + retval.inv_dist_ = 0.0; + retval.factor_ = 0.0; + } + // Else search in the vector + else { + //First element greater than or equal to value + //Start with the second element, so that floor_iter does not go out of range + //Don't access out-of-range, therefore values.end()-1 + auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, value); + + //Find last element smaller than range + auto floor_iter = ceil_iter-1; + + //Find the indices + retval.ind_[0] = floor_iter - values.begin(); + retval.ind_[1] = ceil_iter - values.begin(); + + //Find interpolation ratio + double dist = (*ceil_iter - *floor_iter); + if (std::abs(dist) > 0.0) { + //Possible source for floating point error here if value and floor are large, + //but very close to each other + retval.inv_dist_ = 1.0 / dist; + retval.factor_ = (value-*floor_iter) * retval.inv_dist_; + } + else { + retval.inv_dist_ = 0.0; + retval.factor_ = 0.0; + } + } + + return retval; +} + + + + + + + +/** + * Helper function which interpolates data using the indices etc. given in the inputs. + */ +#ifdef __GNUC__ +#pragma GCC push_options +#pragma GCC optimize ("unroll-loops") +#endif +inline VFPProdProperties::adb_like interpolate( + const VFPProdTable::array_type& array, + const InterpData& flo_i, + const InterpData& thp_i, + const InterpData& wfr_i, + const InterpData& gfr_i, + const InterpData& alq_i) { + + //Values and derivatives in a 5D hypercube + VFPProdProperties::adb_like nn[2][2][2][2][2]; + + + //Pick out nearest neighbors (nn) to our evaluation point + //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. + //The following ladder of for loops will presumably be unrolled by a reasonable compiler. + 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 + const int ti = thp_i.ind_[t]; + const int wi = wfr_i.ind_[w]; + const int gi = gfr_i.ind_[g]; + const int ai = alq_i.ind_[a]; + const int fi = flo_i.ind_[f]; + + //Copy element + nn[t][w][g][a][f].value = array[ti][wi][gi][ai][fi]; + } + } + } + } + } + + //Calculate derivatives + //Note that the derivative of the two end points of a line aligned with the + //"axis of the derivative" are equal + for (int i=0; i<=1; ++i) { + for (int j=0; j<=1; ++j) { + for (int k=0; k<=1; ++k) { + for (int l=0; l<=1; ++l) { + nn[0][i][j][k][l].dthp = (nn[1][i][j][k][l].value - nn[0][i][j][k][l].value) * thp_i.inv_dist_; + nn[i][0][j][k][l].dwfr = (nn[i][1][j][k][l].value - nn[i][0][j][k][l].value) * wfr_i.inv_dist_; + nn[i][j][0][k][l].dgfr = (nn[i][j][1][k][l].value - nn[i][j][0][k][l].value) * gfr_i.inv_dist_; + nn[i][j][k][0][l].dalq = (nn[i][j][k][1][l].value - nn[i][j][k][0][l].value) * alq_i.inv_dist_; + nn[i][j][k][l][0].dflo = (nn[i][j][k][l][1].value - nn[i][j][k][l][0].value) * flo_i.inv_dist_; + + nn[1][i][j][k][l].dthp = nn[0][i][j][k][l].dthp; + nn[i][1][j][k][l].dwfr = nn[i][0][j][k][l].dwfr; + nn[i][j][1][k][l].dgfr = nn[i][j][0][k][l].dgfr; + nn[i][j][k][1][l].dalq = nn[i][j][k][0][l].dalq; + nn[i][j][k][l][1].dflo = nn[i][j][k][l][0].dflo; + } + } + } + } + + double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. + + // Remove dimensions one by one + // 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. + t2 = flo_i.factor_; + t1 = (1.0-t2); + 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] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1]; + } + } + } + } + + t2 = alq_i.factor_; + t1 = (1.0-t2); + 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] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0]; + } + } + } + + t2 = gfr_i.factor_; + t1 = (1.0-t2); + for (int t=0; t<=1; ++t) { + for (int w=0; w<=1; ++w) { + nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0]; + } + } + + t2 = wfr_i.factor_; + t1 = (1.0-t2); + for (int t=0; t<=1; ++t) { + nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0]; + } + + t2 = thp_i.factor_; + t1 = (1.0-t2); + nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0]; + + return nn[0][0][0][0][0]; +} + +#ifdef __GNUC__ +#pragma GCC pop_options //unroll loops +#endif + + + + + +} // namespace detail + + +} // namespace + + + + +#endif /* OPM_AUTODIFF_VFPHELPERS_HPP_ */ diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index bd02e04e5..c17fe1fd4 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -23,57 +23,13 @@ #include #include +#include + namespace Opm { - - - - -VFPProdProperties::VFPProdProperties() { - -} - -VFPProdProperties::VFPProdProperties(const VFPProdTable* table){ - m_tables[table->getTableNum()] = table; -} - - -VFPProdProperties::VFPProdProperties(const std::map& tables) { - init(tables); -} - - - - - -void VFPProdProperties::init(const std::map& prod_tables) { - //Populate production table pointers - for (const auto& table : prod_tables) { - m_tables[table.first] = &table.second; - } -} - -VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, - const Wells& wells, - const ADB& qs, - const ADB& thp, - const ADB& alq) const { - const int np = wells.number_of_phases; - const int nw = wells.number_of_wells; - - //Short-hands for water / oil / gas phases - //TODO enable support for two-phase. - assert(np == 3); - const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); - const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); - const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); - - return bhp(table_id, w, o, g, thp, alq); -} - namespace detail { /** * Returns the type variable for FLO/GFR/WFR @@ -111,7 +67,7 @@ namespace detail { const VFPProdProperties::ADB& liquid, const VFPProdProperties::ADB& vapour, VFPProdTable::FLO_TYPE type) { - return VFPProdProperties::getFlo(aqua, liquid, vapour, type); + return detail::getFlo(aqua, liquid, vapour, type); } template <> @@ -120,7 +76,7 @@ namespace detail { const VFPProdProperties::ADB& liquid, const VFPProdProperties::ADB& vapour, VFPProdTable::WFR_TYPE type) { - return VFPProdProperties::getWFR(aqua, liquid, vapour, type); + return detail::getWFR(aqua, liquid, vapour, type); } template <> @@ -129,7 +85,7 @@ namespace detail { const VFPProdProperties::ADB& liquid, const VFPProdProperties::ADB& vapour, VFPProdTable::GFR_TYPE type) { - return VFPProdProperties::getGFR(aqua, liquid, vapour, type); + return detail::getGFR(aqua, liquid, vapour, type); } /** @@ -200,6 +156,12 @@ namespace detail { return retval; } + + + + /** + * Sets block_pattern to be the "union of x.blockPattern() and block_pattern". + */ void extendBlockPattern(const VFPProdProperties::ADB& x, std::vector& block_pattern) { std::vector x_block_pattern = x.blockPattern(); @@ -219,6 +181,9 @@ namespace detail { } } + /** + * Finds the common block pattern for all inputs + */ std::vector commonBlockPattern( const VFPProdProperties::ADB& x1, const VFPProdProperties::ADB& x2, @@ -236,8 +201,109 @@ namespace detail { return block_pattern; } + /** + * Helper function that finds x for a given value of y for a line + * *NOTE ORDER OF ARGUMENTS* + */ + double findX(const double& x0, + const double& x1, + const double& y0, + const double& y1, + const double& y) { + const double dx = x1 - x0; + const double dy = y1 - y0; + + /** + * y = y0 + (dy / dx) * (x - x0) + * => x = x0 + (y - y0) * (dx / dy) + * + * If dy is zero, use x1 as the value. + */ + + double x = 0.0; + + if (dy != 0.0) { + x = x0 + (y-y0) * (dx/dy); + } + else { + x = x1; + } + + return x; + } } //Namespace + + + + + + + + + + + + + + + + + + + + +VFPProdProperties::VFPProdProperties() { + +} + + + +VFPProdProperties::VFPProdProperties(const VFPProdTable* table){ + m_tables[table->getTableNum()] = table; +} + + + + +VFPProdProperties::VFPProdProperties(const std::map& tables) { + init(tables); +} + + + +void VFPProdProperties::init(const std::map& prod_tables) { + //Populate production table pointers + for (const auto& table : prod_tables) { + m_tables[table.first] = &table.second; + } +} + + + +VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, + const Wells& wells, + const ADB& qs, + const ADB& thp, + const ADB& alq) const { + const int np = wells.number_of_phases; + const int nw = wells.number_of_wells; + + //Short-hands for water / oil / gas phases + //TODO enable support for two-phase. + assert(np == 3); + const ADB& w = subset(qs, Span(nw, 1, BlackoilPhases::Aqua*nw)); + const ADB& o = subset(qs, Span(nw, 1, BlackoilPhases::Liquid*nw)); + const ADB& g = subset(qs, Span(nw, 1, BlackoilPhases::Vapour*nw)); + + return bhp(table_id, w, o, g, thp, alq); +} + + + + + + VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, const ADB& aqua, const ADB& liquid, @@ -281,13 +347,13 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, const VFPProdTable* table = well_tables[i]; if (table != NULL) { //First, find the values to interpolate between - auto flo_i = find_interp_data(flo.value()[i], table->getFloAxis()); - auto thp_i = find_interp_data(thp.value()[i], table->getTHPAxis()); - auto wfr_i = find_interp_data(wfr.value()[i], table->getWFRAxis()); - auto gfr_i = find_interp_data(gfr.value()[i], table->getGFRAxis()); - auto alq_i = find_interp_data(alq.value()[i], table->getALQAxis()); + 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()); - adb_like bhp_val = interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + adb_like bhp_val = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); value[i] = bhp_val.value; dthp[i] = bhp_val.dthp; @@ -337,6 +403,9 @@ VFPProdProperties::ADB VFPProdProperties::bhp(const std::vector& table_id, return retval; } + + + VFPProdProperties::adb_like VFPProdProperties::bhp(int table_id, const double& aqua, const double& liquid, @@ -346,23 +415,26 @@ VFPProdProperties::adb_like VFPProdProperties::bhp(int table_id, const VFPProdTable* table = getProdTable(table_id); //Find interpolation variables - double flo = getFlo(aqua, liquid, vapour, table->getFloType()); - double wfr = getWFR(aqua, liquid, vapour, table->getWFRType()); - double gfr = getGFR(aqua, liquid, vapour, table->getGFRType()); + double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); + double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); //First, find the values to interpolate between - auto flo_i = find_interp_data(flo, table->getFloAxis()); - auto thp_i = find_interp_data(thp, table->getTHPAxis()); - auto wfr_i = find_interp_data(wfr, table->getWFRAxis()); - auto gfr_i = find_interp_data(gfr, table->getGFRAxis()); - auto alq_i = find_interp_data(alq, table->getALQAxis()); + 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()); //Then perform the interpolation itself - adb_like retval = interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); + adb_like retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i); return retval; } + + + double VFPProdProperties::thp(int table_id, const double& aqua, const double& liquid, @@ -375,9 +447,9 @@ double VFPProdProperties::thp(int table_id, double thp = -1e100; //Find interpolation variables - double flo = getFlo(aqua, liquid, vapour, table->getFloType()); - double wfr = getWFR(aqua, liquid, vapour, table->getWFRType()); - double gfr = getGFR(aqua, liquid, vapour, table->getGFRType()); + double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType()); + double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType()); + double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType()); /** * Get THP axis, assume that it is sorted @@ -391,14 +463,14 @@ double VFPProdProperties::thp(int table_id, * by interpolating for every value of thp. This might be somewhat * expensive, but let us assome that nthp is small */ - auto flo_i = find_interp_data(flo, table->getFloAxis()); - auto wfr_i = find_interp_data(wfr, table->getWFRAxis()); - auto gfr_i = find_interp_data(gfr, table->getGFRAxis()); - auto alq_i = find_interp_data(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 bhp_array[nthp-1]) { @@ -424,7 +496,7 @@ double VFPProdProperties::thp(int table_id, const double& x1 = thp_array[nthp-1]; const double& y0 = bhp_array[nthp-2]; const double& y1 = bhp_array[nthp-1]; - thp = find_x(x0, x1, y0, y1, bhp); + thp = detail::findX(x0, x1, y0, y1, bhp); } //Target bhp within table ranges, interpolate else { @@ -452,7 +524,7 @@ double VFPProdProperties::thp(int table_id, const double& x1 = thp_array[i+1]; const double& y0 = bhp_array[i ]; const double& y1 = bhp_array[i+1]; - thp = find_x(x0, x1, y0, y1, bhp); + thp = detail::findX(x0, x1, y0, y1, bhp); } } //bhp_array not sorted, raw search. @@ -476,7 +548,7 @@ double VFPProdProperties::thp(int table_id, const double& x1 = thp_array[i+1]; const double& y0 = bhp_array[i ]; const double& y1 = bhp_array[i+1]; - thp = find_x(x0, x1, y0, y1, bhp); + thp = detail::findX(x0, x1, y0, y1, bhp); } else if (bhp <= bhp_array[0]) { //TODO: LOG extrapolation @@ -484,7 +556,7 @@ double VFPProdProperties::thp(int table_id, const double& x1 = thp_array[1]; const double& y0 = bhp_array[0]; const double& y1 = bhp_array[1]; - thp = find_x(x0, x1, y0, y1, bhp); + thp = detail::findX(x0, x1, y0, y1, bhp); } //Target bhp greater than all values in array, extrapolate else if (bhp > bhp_array[nthp-1]) { @@ -493,7 +565,7 @@ double VFPProdProperties::thp(int table_id, const double& x1 = thp_array[nthp-1]; const double& y0 = bhp_array[nthp-2]; const double& y1 = bhp_array[nthp-1]; - thp = find_x(x0, x1, y0, y1, bhp); + thp = detail::findX(x0, x1, y0, y1, bhp); } else { OPM_THROW(std::logic_error, "Programmer error: Unable to find THP in THP array"); @@ -505,6 +577,9 @@ double VFPProdProperties::thp(int table_id, + + + const VFPProdTable* VFPProdProperties::getProdTable(int table_id) const { auto entry = m_tables.find(table_id); if (entry == m_tables.end()) { @@ -515,193 +590,10 @@ const VFPProdTable* VFPProdProperties::getProdTable(int table_id) const { } } -VFPProdProperties::InterpData VFPProdProperties::find_interp_data(const double& value, const std::vector& values) { - InterpData retval; - - //If we only have one value in our vector, return that - if (values.size() == 1) { - retval.ind_[0] = 0; - retval.ind_[1] = 0; - retval.inv_dist_ = 0.0; - retval.factor_ = 0.0; - } - // Else search in the vector - else { - //First element greater than or equal to value - //Start with the second element, so that floor_iter does not go out of range - //Don't access out-of-range, therefore values.end()-1 - auto ceil_iter = std::lower_bound(values.begin()+1, values.end()-1, value); - - //Find last element smaller than range - auto floor_iter = ceil_iter-1; - - //Find the indices - retval.ind_[0] = floor_iter - values.begin(); - retval.ind_[1] = ceil_iter - values.begin(); - - //Find interpolation ratio - double dist = (*ceil_iter - *floor_iter); - if (std::abs(dist) > 0.0) { - //Possible source for floating point error here if value and floor are large, - //but very close to each other - retval.inv_dist_ = 1.0 / dist; - retval.factor_ = (value-*floor_iter) * retval.inv_dist_; - } - else { - retval.inv_dist_ = 0.0; - retval.factor_ = 0.0; - } - } - - return retval; -} -#ifdef __GNUC__ -#pragma GCC push_options -#pragma GCC optimize ("unroll-loops") -#endif - -VFPProdProperties::adb_like VFPProdProperties::interpolate( - const VFPProdTable::array_type& array, - const InterpData& flo_i, - const InterpData& thp_i, - const InterpData& wfr_i, - const InterpData& gfr_i, - const InterpData& alq_i) { - - //Values and derivatives in a 5D hypercube - adb_like nn[2][2][2][2][2]; - //Pick out nearest neighbors (nn) to our evaluation point - //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. - //The following ladder of for loops will presumably be unrolled by a reasonable compiler. - 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 - const int ti = thp_i.ind_[t]; - const int wi = wfr_i.ind_[w]; - const int gi = gfr_i.ind_[g]; - const int ai = alq_i.ind_[a]; - const int fi = flo_i.ind_[f]; - - //Copy element - nn[t][w][g][a][f].value = array[ti][wi][gi][ai][fi]; - } - } - } - } - } - - //Calculate derivatives - //Note that the derivative of the two end points of a line aligned with the - //"axis of the derivative" are equal - for (int i=0; i<=1; ++i) { - for (int j=0; j<=1; ++j) { - for (int k=0; k<=1; ++k) { - for (int l=0; l<=1; ++l) { - nn[0][i][j][k][l].dthp = (nn[1][i][j][k][l].value - nn[0][i][j][k][l].value) * thp_i.inv_dist_; - nn[i][0][j][k][l].dwfr = (nn[i][1][j][k][l].value - nn[i][0][j][k][l].value) * wfr_i.inv_dist_; - nn[i][j][0][k][l].dgfr = (nn[i][j][1][k][l].value - nn[i][j][0][k][l].value) * gfr_i.inv_dist_; - nn[i][j][k][0][l].dalq = (nn[i][j][k][1][l].value - nn[i][j][k][0][l].value) * alq_i.inv_dist_; - nn[i][j][k][l][0].dflo = (nn[i][j][k][l][1].value - nn[i][j][k][l][0].value) * flo_i.inv_dist_; - - nn[1][i][j][k][l].dthp = nn[0][i][j][k][l].dthp; - nn[i][1][j][k][l].dwfr = nn[i][0][j][k][l].dwfr; - nn[i][j][1][k][l].dgfr = nn[i][j][0][k][l].dgfr; - nn[i][j][k][1][l].dalq = nn[i][j][k][0][l].dalq; - nn[i][j][k][l][1].dflo = nn[i][j][k][l][0].dflo; - } - } - } - } - - double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. - - // Remove dimensions one by one - // 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. - t2 = flo_i.factor_; - t1 = (1.0-t2); - 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] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1]; - } - } - } - } - - t2 = alq_i.factor_; - t1 = (1.0-t2); - 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] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0]; - } - } - } - - t2 = gfr_i.factor_; - t1 = (1.0-t2); - for (int t=0; t<=1; ++t) { - for (int w=0; w<=1; ++w) { - nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0]; - } - } - - t2 = wfr_i.factor_; - t1 = (1.0-t2); - for (int t=0; t<=1; ++t) { - nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0]; - } - - t2 = thp_i.factor_; - t1 = (1.0-t2); - nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0]; - - return nn[0][0][0][0][0]; -} - -#ifdef __GNUC__ -#pragma GCC pop_options //unroll loops -#endif - - -double VFPProdProperties::find_x(const double& x0, - const double& x1, - const double& y0, - const double& y1, - const double& y) { - const double dx = x1 - x0; - const double dy = y1 - y0; - - /** - * y = y0 + (dy / dx) * (x - x0) - * => x = x0 + (y - y0) * (dx / dy) - * - * If dy is zero, use x1 as the value. - */ - - double x = 0.0; - - if (dy != 0.0) { - x = x0 + (y-y0) * (dx/dy); - } - else { - x = x1; - } - - return x; -} diff --git a/opm/autodiff/VFPProdProperties.hpp b/opm/autodiff/VFPProdProperties.hpp index d62723687..d93d3c239 100644 --- a/opm/autodiff/VFPProdProperties.hpp +++ b/opm/autodiff/VFPProdProperties.hpp @@ -152,117 +152,11 @@ public: const double& bhp, const double& alq) const; - /** - * Computes the flo parameter according to the flo_type_ - * @return Production rate of oil, gas or liquid. - */ - template - static T getFlo(const T& aqua, const T& liquid, const T& vapour, - const VFPProdTable::FLO_TYPE& type) { - switch (type) { - case VFPProdTable::FLO_OIL: - //Oil = liquid phase - return liquid; - case VFPProdTable::FLO_LIQ: - //Liquid = aqua + liquid phases - return aqua + liquid; - case VFPProdTable::FLO_GAS: - //Gas = vapor phase - return vapour; - case VFPProdTable::FLO_INVALID: //Intentional fall-through - default: - OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'"); - } - } - - - /** - * Computes the wfr parameter according to the wfr_type_ - * @return Production rate of oil, gas or liquid. - */ - template - static T getWFR(const T& aqua, const T& liquid, const T& vapour, - const VFPProdTable::WFR_TYPE& type) { - switch(type) { - case VFPProdTable::WFR_WOR: { - //Water-oil ratio = water / oil - T wor = aqua / liquid; - return zeroIfNan(wor); - } - case VFPProdTable::WFR_WCT: - //Water cut = water / (water + oil) - return zeroIfNan(aqua / (aqua + liquid)); - case VFPProdTable::WFR_WGR: - //Water-gas ratio = water / gas - return zeroIfNan(aqua / vapour); - case VFPProdTable::WFR_INVALID: //Intentional fall-through - default: - OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'"); - } - } - - /** - * Computes the gfr parameter according to the gfr_type_ - * @return Production rate of oil, gas or liquid. - */ - template - static T getGFR(const T& aqua, const T& liquid, const T& vapour, - const VFPProdTable::GFR_TYPE& type) { - switch(type) { - case VFPProdTable::GFR_GOR: - // Gas-oil ratio = gas / oil - return zeroIfNan(vapour / liquid); - case VFPProdTable::GFR_GLR: - // Gas-liquid ratio = gas / (oil + water) - return zeroIfNan(vapour / (liquid + aqua)); - case VFPProdTable::GFR_OGR: - // Oil-gas ratio = oil / gas - return zeroIfNan(liquid / vapour); - case VFPProdTable::GFR_INVALID: //Intentional fall-through - default: - OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'"); - } - } - private: // Map which connects the table number with the table itself std::map m_tables; - /** - * Helper struct for linear interpolation - */ - struct InterpData { - InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {} - int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value] - double inv_dist_; // 1 / distance between the two end points of the segment. Used to calculate derivatives and uses 1.0 / 0.0 = 0.0 as a convention - double factor_; // Interpolation factor - }; - - /** - * Helper function to find indices etc. for linear interpolation - */ - static InterpData find_interp_data(const double& value, const std::vector& values); - - /** - * Helper function which interpolates data using the indices etc. given in the inputs. - */ - static adb_like interpolate(const VFPProdTable::array_type& array, - const InterpData& flo_i, - const InterpData& thp_i, - const InterpData& wfr_i, - const InterpData& gfr_i, - const InterpData& alq_i); - - /** - * Helper function that finds x for a given value of y for a line - * *NOTE ORDER OF ARGUMENTS* - */ - static double find_x(const double& x0, - const double& x1, - const double& y0, - const double& y1, - const double& y); /** @@ -275,19 +169,6 @@ private: */ const VFPProdTable* getProdTable(int table_id) const; - static inline double zeroIfNan(const double& value) { - return (std::isnan(value)) ? 0.0 : value; - } - - static inline ADB zeroIfNan(const ADB& values) { - Selector not_nan_selector(values.value(), Selector::NotNaN); - - const ADB::V z = ADB::V::Zero(values.size()); - const ADB zero = ADB::constant(z, values.blockPattern()); - - ADB retval = not_nan_selector.select(values, zero); - return retval; - } }; diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index b293e26ba..c77b81e9e 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -42,11 +42,182 @@ #include #include +#include const double max_d_tol = 1.0e-10; const double sad_tol = 1.0e-8; + + + + + + + + + +struct ConversionFixture { + typedef Opm::VFPProdProperties::ADB ADB; + + ConversionFixture() : + num_wells(5), + aqua(ADB::null()), + liquid(ADB::null()), + vapour(ADB::null()) + { + ADB::V aqua_v(num_wells); + ADB::V liquid_v(num_wells); + ADB::V vapour_v(num_wells); + + for (int i=0; i ref_flo_oil(num_wells); + std::vector ref_flo_liq(num_wells); + std::vector ref_flo_gas(num_wells); + for (int i=0; i ref_wfr_wor(num_wells); + std::vector ref_wfr_wct(num_wells); + std::vector ref_wfr_wgr(num_wells); + for (int i=0; i ref_gfr_gor(num_wells); + std::vector ref_gfr_glr(num_wells); + std::vector ref_gfr_ogr(num_wells); + for (int i=0; i(n); //Find values that should be in table - double flo = properties->getFlo(aqua, liquid, vapour, table.getFloType()); - double wfr = properties->getWFR(aqua, liquid, vapour, table.getWFRType()); - double gfr = properties->getGFR(aqua, liquid, vapour, table.getGFRType()); + double flo = Opm::detail::getFlo(aqua, liquid, vapour, table.getFloType()); + double wfr = Opm::detail::getWFR(aqua, liquid, vapour, table.getWFRType()); + double gfr = Opm::detail::getGFR(aqua, liquid, vapour, table.getGFRType()); //Calculate reference adb_like reference; @@ -437,9 +608,9 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlane) const double liquid = m / static_cast(n); //Find values that should be in table - double v = properties->getFlo(aqua, liquid, vapour, table.getFloType()); - double y = properties->getWFR(aqua, liquid, vapour, table.getWFRType()); - double z = properties->getGFR(aqua, liquid, vapour, table.getGFRType()); + 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; reference_sum += reference; @@ -524,9 +695,9 @@ BOOST_AUTO_TEST_CASE(ExtrapolatePlaneADB) double reference = 0.0; for (int w=0; w < num_wells; ++w) { //Find values that should be in table - double v = properties->getFlo(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getFloType()); - double y = properties->getWFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getWFRType()); - double z = properties->getGFR(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getGFRType()); + double v = Opm::detail::getFlo(aqua*(w+1), liquid*(w+1), vapour*(w+1), table.getFloType()); + 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; value = bhp_val[w]; @@ -668,155 +839,6 @@ BOOST_AUTO_TEST_SUITE_END() // Trivial tests -struct ConversionFixture { - typedef Opm::VFPProdProperties::ADB ADB; - - ConversionFixture() : - num_wells(5), - aqua(ADB::null()), - liquid(ADB::null()), - vapour(ADB::null()) - { - ADB::V aqua_v(num_wells); - ADB::V liquid_v(num_wells); - ADB::V vapour_v(num_wells); - - for (int i=0; i ref_flo_oil(num_wells); - std::vector ref_flo_liq(num_wells); - std::vector ref_flo_gas(num_wells); - for (int i=0; i ref_wfr_wor(num_wells); - std::vector ref_wfr_wct(num_wells); - std::vector ref_wfr_wgr(num_wells); - for (int i=0; i ref_gfr_gor(num_wells); - std::vector ref_gfr_glr(num_wells); - std::vector ref_gfr_ogr(num_wells); - for (int i=0; i