diff --git a/examples/printvfp.cpp b/examples/printvfp.cpp index a58fd6b3b..698b01478 100644 --- a/examples/printvfp.cpp +++ b/examples/printvfp.cpp @@ -83,13 +83,13 @@ double computeBhp(const VFPProdTable& table, // First, find the values to interpolate between. // Assuming positive flo here! assert(flo > 0.0); - auto flo_i = detail::findInterpData(flo, table.getFloAxis()); - auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); // assume constant - auto wfr_i = detail::findInterpData(wfr, table.getWFRAxis()); - auto gfr_i = detail::findInterpData(gfr, table.getGFRAxis()); - auto alq_i = detail::findInterpData(alq, table.getALQAxis()); //assume constant + auto flo_i = VFPHelpers::findInterpData(flo, table.getFloAxis()); + auto thp_i = VFPHelpers::findInterpData(thp, table.getTHPAxis()); // assume constant + auto wfr_i = VFPHelpers::findInterpData(wfr, table.getWFRAxis()); + auto gfr_i = VFPHelpers::findInterpData(gfr, table.getGFRAxis()); + auto alq_i = VFPHelpers::findInterpData(alq, table.getALQAxis()); //assume constant - return detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; + return VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; } diff --git a/opm/simulators/wells/VFPHelpers.cpp b/opm/simulators/wells/VFPHelpers.cpp index c7c847263..ccb859b33 100644 --- a/opm/simulators/wells/VFPHelpers.cpp +++ b/opm/simulators/wells/VFPHelpers.cpp @@ -38,14 +38,15 @@ namespace { * 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) +template +Scalar findX(const Scalar x0, + const Scalar x1, + const Scalar y0, + const Scalar y1, + const Scalar y) { - const double dx = x1 - x0; - const double dy = y1 - y0; + const Scalar dx = x1 - x0; + const Scalar dy = y1 - y0; /** * y = y0 + (dy / dx) * (x - x0) @@ -54,7 +55,7 @@ double findX(const double x0, * If dy is zero, use x1 as the value. */ - double x = 0.0; + Scalar x = 0.0; if (dy != 0.0) { x = x0 + (y-y0) * (dx/dy); @@ -77,17 +78,18 @@ static T chopNegativeValues(const T& value) { } namespace Opm { -namespace detail { -InterpData findInterpData(const double value_in, const std::vector& values) +template +detail::InterpData VFPHelpers::findInterpData(const Scalar value_in, + const std::vector& values) { - InterpData retval; + detail::InterpData retval; const int nvalues = values.size(); // chopping the value to be zero, which means we do not // extrapolate the table towards nagative ranges - const double value = value_in < 0.? 0. : value_in; + const Scalar value = value_in < 0.? 0. : value_in; //If we only have one value in our vector, return that if (values.size() == 1) { @@ -119,8 +121,8 @@ InterpData findInterpData(const double value_in, const std::vector& valu } } - const double start = values[retval.ind_[0]]; - const double end = values[retval.ind_[1]]; + const Scalar start = values[retval.ind_[0]]; + const Scalar end = values[retval.ind_[1]]; //Find interpolation ratio if (end > start) { @@ -138,50 +140,17 @@ InterpData findInterpData(const double value_in, const std::vector& valu return retval; } -VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs) -{ - lhs.value += rhs.value; - lhs.dthp += rhs.dthp; - lhs.dwfr += rhs.dwfr; - lhs.dgfr += rhs.dgfr; - lhs.dalq += rhs.dalq; - lhs.dflo += rhs.dflo; - return lhs; -} - -VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs) -{ - lhs.value -= rhs.value; - lhs.dthp -= rhs.dthp; - lhs.dwfr -= rhs.dwfr; - lhs.dgfr -= rhs.dgfr; - lhs.dalq -= rhs.dalq; - lhs.dflo -= rhs.dflo; - return lhs; -} - -VFPEvaluation operator*(double lhs, const VFPEvaluation& rhs) -{ - VFPEvaluation retval; - retval.value = rhs.value * lhs; - retval.dthp = rhs.dthp * lhs; - retval.dwfr = rhs.dwfr * lhs; - retval.dgfr = rhs.dgfr * lhs; - retval.dalq = rhs.dalq * lhs; - retval.dflo = rhs.dflo * lhs; - return retval; -} - -VFPEvaluation interpolate(const VFPProdTable& table, - const InterpData& flo_i, - const InterpData& thp_i, - const InterpData& wfr_i, - const InterpData& gfr_i, - const InterpData& alq_i) +template +detail::VFPEvaluation VFPHelpers:: +interpolate(const VFPProdTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i, + const detail::InterpData& wfr_i, + const detail::InterpData& gfr_i, + const detail::InterpData& alq_i) { //Values and derivatives in a 5D hypercube - VFPEvaluation nn[2][2][2][2][2]; - + detail::VFPEvaluation 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 @@ -231,7 +200,7 @@ VFPEvaluation interpolate(const VFPProdTable& table, } } - double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. + Scalar 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 @@ -280,13 +249,14 @@ VFPEvaluation interpolate(const VFPProdTable& table, return nn[0][0][0][0][0]; } -VFPEvaluation interpolate(const VFPInjTable& table, - const InterpData& flo_i, - const InterpData& thp_i) +template +detail::VFPEvaluation VFPHelpers:: +interpolate(const VFPInjTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i) { //Values and derivatives in a 2D plane - VFPEvaluation nn[2][2]; - + detail::VFPEvaluation nn[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. @@ -318,7 +288,7 @@ VFPEvaluation interpolate(const VFPInjTable& table, nn[i][1].dflo = nn[i][0].dflo; } - double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. + Scalar t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t. // Go from 2D to 1D t2 = flo_i.factor_; @@ -334,20 +304,22 @@ VFPEvaluation interpolate(const VFPInjTable& table, return nn[0][0]; } -VFPEvaluation bhp(const VFPProdTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp, - const double alq, - const double explicit_wfr, - const double explicit_gfr, - const bool use_vfpexplicit) +template +detail::VFPEvaluation VFPHelpers:: +bhp(const VFPProdTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_vfpexplicit) { //Find interpolation variables - double flo = detail::getFlo(table, aqua, liquid, vapour); - double wfr = detail::getWFR(table, aqua, liquid, vapour); - double gfr = detail::getGFR(table, aqua, liquid, vapour); + Scalar flo = detail::getFlo(table, aqua, liquid, vapour); + Scalar wfr = detail::getWFR(table, aqua, liquid, vapour); + Scalar gfr = detail::getGFR(table, aqua, liquid, vapour); if (use_vfpexplicit || -flo < table.getFloAxis().front()) { wfr = explicit_wfr; gfr = explicit_gfr; @@ -355,43 +327,47 @@ VFPEvaluation bhp(const VFPProdTable& table, //First, find the values to interpolate between //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()); + auto flo_i = findInterpData(-flo, table.getFloAxis()); + auto thp_i = findInterpData( thp, table.getTHPAxis()); + auto wfr_i = findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = findInterpData( gfr, table.getGFRAxis()); + auto alq_i = findInterpData( alq, table.getALQAxis()); - detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + detail::VFPEvaluation retval = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); return retval; } -VFPEvaluation bhp(const VFPInjTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp) +template +detail::VFPEvaluation VFPHelpers:: +bhp(const VFPInjTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp) { //Find interpolation variables - double flo = detail::getFlo(table, aqua, liquid, vapour); + Scalar flo = detail::getFlo(table, aqua, liquid, vapour); //First, find the values to interpolate between - auto flo_i = detail::findInterpData(flo, table.getFloAxis()); - auto thp_i = detail::findInterpData(thp, table.getTHPAxis()); + auto flo_i = findInterpData(flo, table.getFloAxis()); + auto thp_i = findInterpData(thp, table.getTHPAxis()); //Then perform the interpolation itself - detail::VFPEvaluation retval = detail::interpolate(table, flo_i, thp_i); + detail::VFPEvaluation retval = interpolate(table, flo_i, thp_i); return retval; } -double findTHP(const std::vector& bhp_array, - const std::vector& thp_array, - double bhp) +template +Scalar VFPHelpers:: +findTHP(const std::vector& bhp_array, + const std::vector& thp_array, + Scalar bhp) { int nthp = thp_array.size(); - double thp = -1e100; + Scalar thp = -1e100; //Check that our thp axis is sorted assert(std::is_sorted(thp_array.begin(), thp_array.end())); @@ -406,19 +382,19 @@ double findTHP(const std::vector& bhp_array, //Target bhp less than all values in array, extrapolate if (bhp <= bhp_array[0]) { //TODO: LOG extrapolation - const double& x0 = thp_array[0]; - const double& x1 = thp_array[1]; - const double& y0 = bhp_array[0]; - const double& y1 = bhp_array[1]; + const Scalar& x0 = thp_array[0]; + const Scalar& x1 = thp_array[1]; + const Scalar& y0 = bhp_array[0]; + const Scalar& y1 = bhp_array[1]; thp = findX(x0, x1, y0, y1, bhp); } //Target bhp greater than all values in array, extrapolate else if (bhp > bhp_array[nthp-1]) { //TODO: LOG extrapolation - const double& x0 = thp_array[nthp-2]; - const double& x1 = thp_array[nthp-1]; - const double& y0 = bhp_array[nthp-2]; - const double& y1 = bhp_array[nthp-1]; + const Scalar& x0 = thp_array[nthp-2]; + const Scalar& x1 = thp_array[nthp-1]; + const Scalar& y0 = bhp_array[nthp-2]; + const Scalar& y1 = bhp_array[nthp-1]; thp = findX(x0, x1, y0, y1, bhp); } //Target bhp within table ranges, interpolate @@ -432,8 +408,8 @@ double findTHP(const std::vector& bhp_array, int i=0; bool found = false; for (; i& bhp_array, assert(found == true); static_cast(found); //Silence compiler warning - const double& x0 = thp_array[i ]; - const double& x1 = thp_array[i+1]; - const double& y0 = bhp_array[i ]; - const double& y1 = bhp_array[i+1]; + const Scalar& x0 = thp_array[i ]; + const Scalar& x1 = thp_array[i+1]; + const Scalar& y0 = bhp_array[i ]; + const Scalar& y1 = bhp_array[i+1]; thp = findX(x0, x1, y0, y1, bhp); } } @@ -459,8 +435,8 @@ double findTHP(const std::vector& bhp_array, int i=0; bool found = false; for (; i& bhp_array, } } if (found) { - const double& x0 = thp_array[i ]; - const double& x1 = thp_array[i+1]; - const double& y0 = bhp_array[i ]; - const double& y1 = bhp_array[i+1]; + const Scalar& x0 = thp_array[i ]; + const Scalar& x1 = thp_array[i+1]; + const Scalar& y0 = bhp_array[i ]; + const Scalar& y1 = bhp_array[i+1]; thp = findX(x0, x1, y0, y1, bhp); } else if (bhp <= bhp_array[0]) { //TODO: LOG extrapolation - const double& x0 = thp_array[0]; - const double& x1 = thp_array[1]; - const double& y0 = bhp_array[0]; - const double& y1 = bhp_array[1]; + const Scalar& x0 = thp_array[0]; + const Scalar& x1 = thp_array[1]; + const Scalar& y0 = bhp_array[0]; + const Scalar& y1 = bhp_array[1]; thp = findX(x0, x1, y0, y1, bhp); } //Target bhp greater than all values in array, extrapolate else if (bhp > bhp_array[nthp-1]) { //TODO: LOG extrapolation - const double& x0 = thp_array[nthp-2]; - const double& x1 = thp_array[nthp-1]; - const double& y0 = bhp_array[nthp-2]; - const double& y1 = bhp_array[nthp-1]; + const Scalar& x0 = thp_array[nthp-2]; + const Scalar& x1 = thp_array[nthp-1]; + const Scalar& y0 = bhp_array[nthp-2]; + const Scalar& y1 = bhp_array[nthp-1]; thp = findX(x0, x1, y0, y1, bhp); } else { @@ -499,86 +475,88 @@ double findTHP(const std::vector& bhp_array, return thp; } -std::pair +template +std::pair VFPHelpers:: getMinimumBHPCoordinate(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq) -{ + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq) +{ // Given fixed thp, wfr, gfr and alq, this function finds the minimum bhp and returns - // the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the - // shape of the function bhp(flo), so all points in the flo-axis is checked. - double flo_at_bhp_min = 0.0; // start by checking flo=0 - auto flo_i = detail::findInterpData(flo_at_bhp_min, 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()); + // the corresponding pair (-flo_at_bhp_min, bhp_min). No assumption is taken on the + // shape of the function bhp(flo), so all points in the flo-axis is checked. + Scalar flo_at_bhp_min = 0.0; // start by checking flo=0 + auto flo_i = findInterpData(flo_at_bhp_min, table.getFloAxis()); + auto thp_i = findInterpData( thp, table.getTHPAxis()); + auto wfr_i = findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = findInterpData( gfr, table.getGFRAxis()); + auto alq_i = findInterpData( alq, table.getALQAxis()); - detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); - double bhp_min = bhp_i.value; + detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + Scalar bhp_min = bhp_i.value; const std::vector& flos = table.getFloAxis(); for (size_t i = 0; i < flos.size(); ++i) { - flo_i = detail::findInterpData(flos[i], flos); - bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + flo_i = findInterpData(flos[i], flos); + bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); if (bhp_i.value < bhp_min){ bhp_min = bhp_i.value; flo_at_bhp_min = flos[i]; } } // return negative flo - return std::make_pair(-flo_at_bhp_min, bhp_min); -} + return std::make_pair(-flo_at_bhp_min, bhp_min); +} -std::optional> +template +std::optional> VFPHelpers:: intersectWithIPR(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq, - const double ipr_a, - const double ipr_b, - const std::function& adjust_bhp) -{ - // Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection - // between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq, + const Scalar ipr_a, + const Scalar ipr_b, + const std::function& adjust_bhp) +{ + // Given fixed thp, wfr, gfr and alq, this function finds a stable (-flo, bhp)-intersection + // between the ipr-line and bhp(flo) if such an intersection exists. For multiple stable // intersections, the one corresponding the largest flo is returned. // The adjust_bhp-function is used to adjust the vfp-table bhp-values to actual bhp-values due - // vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments. + // vfp/well ref-depth differences and/or WVFPDP-related pressure adjustments. // NOTE: ipr-line is q=b*bhp - a! // ipr is given for negative flo, so - // flo = -b*bhp + a, i.e., bhp = -(flo-a)/b - 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()); + // flo = -b*bhp + a, i.e., bhp = -(flo-a)/b + auto thp_i = findInterpData( thp, table.getTHPAxis()); + auto wfr_i = findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = findInterpData( gfr, table.getGFRAxis()); + auto alq_i = findInterpData( alq, table.getALQAxis()); if (ipr_b == 0.0) { // this shouldn't happen, but deal with it to be safe - auto flo_i = detail::findInterpData(ipr_a, table.getFloAxis()); - detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + auto flo_i = findInterpData(ipr_a, table.getFloAxis()); + detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); return std::make_pair(-ipr_a, adjust_bhp(bhp_i.value)); } // find largest flo (flo_x) for which y = bhp(flo) + (flo-a)/b = 0 and dy/dflo > 0 - double flo_x = -1.0; - double flo0, flo1; - double y0, y1; + Scalar flo_x = -1.0; + Scalar flo0, flo1; + Scalar y0, y1; flo0 = 0.0; // start by checking flo=0 - auto flo_i = detail::findInterpData(flo0, table.getFloAxis()); - detail::VFPEvaluation bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + auto flo_i = findInterpData(flo0, table.getFloAxis()); + detail::VFPEvaluation bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); y0 = adjust_bhp(bhp_i.value) - ipr_a/ipr_b; // +0.0/ipr_b - std::vector flos = table.getFloAxis(); + const std::vector& flos = table.getFloAxis(); for (size_t i = 0; i < flos.size(); ++i) { flo1 = flos[i]; - flo_i = detail::findInterpData(flo1, table.getFloAxis()); - bhp_i = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + flo_i = findInterpData(flo1, flos); + bhp_i = interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); y1 = adjust_bhp(bhp_i.value) + (flo1 - ipr_a)/ipr_b; if (y0 < 0 && y1 >= 0){ // crossing with positive slope - double w = -y0/(y1-y0); + Scalar w = -y0/(y1-y0); w = std::clamp(w, 0.0, 1.0); // just to be safe (if y0~y1~0) flo_x = flo0 + w*(flo1 - flo0); } @@ -591,7 +569,46 @@ intersectWithIPR(const VFPProdTable& table, } else { return std::nullopt; } -} +} + +namespace detail { + +template +VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs) +{ + lhs.value += rhs.value; + lhs.dthp += rhs.dthp; + lhs.dwfr += rhs.dwfr; + lhs.dgfr += rhs.dgfr; + lhs.dalq += rhs.dalq; + lhs.dflo += rhs.dflo; + return lhs; +} + +template +VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs) +{ + lhs.value -= rhs.value; + lhs.dthp -= rhs.dthp; + lhs.dwfr -= rhs.dwfr; + lhs.dgfr -= rhs.dgfr; + lhs.dalq -= rhs.dalq; + lhs.dflo -= rhs.dflo; + return lhs; +} + +template +VFPEvaluation operator*(Scalar lhs, const VFPEvaluation& rhs) +{ + VFPEvaluation retval; + retval.value = rhs.value * lhs; + retval.dthp = rhs.dthp * lhs; + retval.dwfr = rhs.dwfr * lhs; + retval.dgfr = rhs.dgfr * lhs; + retval.dalq = rhs.dalq * lhs; + retval.dflo = rhs.dflo * lhs; + return retval; +} template T getFlo(const VFPProdTable& table, @@ -753,4 +770,7 @@ INSTANCE(DenseAd::Evaluation) INSTANCE(DenseAd::Evaluation) } // namespace detail + +template class VFPHelpers; + } // namespace Opm diff --git a/opm/simulators/wells/VFPHelpers.hpp b/opm/simulators/wells/VFPHelpers.hpp index df678da98..761035ee5 100644 --- a/opm/simulators/wells/VFPHelpers.hpp +++ b/opm/simulators/wells/VFPHelpers.hpp @@ -21,7 +21,6 @@ #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_ #define OPM_AUTODIFF_VFPHELPERS_HPP_ -#include #include #include #include @@ -37,6 +36,40 @@ class VFPProdTable; namespace detail { +/** + * An "ADB-like" structure with a single value and a set of derivatives + */ +template +struct VFPEvaluation +{ + VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {}; + Scalar value; + Scalar dthp; + Scalar dwfr; + Scalar dgfr; + Scalar dalq; + Scalar dflo; +}; + +template +VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs); +template +VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs); +template +VFPEvaluation operator*(Scalar lhs, const VFPEvaluation& rhs); + +/** + * Helper struct for linear interpolation + */ +template +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] + Scalar 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 + Scalar factor_; // Interpolation factor +}; + /** * Computes the flo parameter according to the flo_type_ * for production tables @@ -79,76 +112,6 @@ T getGFR(const VFPProdTable& table, const T& liquid, const T& vapour); -/** - * 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 and extrapolation - * @param value_in Value to find in values - * @param values Sorted list of values to search for value in. - * @return Data required to find the interpolated value - */ -InterpData findInterpData(const double value_in, const std::vector& values); - -/** - * An "ADB-like" structure with a single value and a set of derivatives - */ -struct VFPEvaluation { - VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {}; - double value; - double dthp; - double dwfr; - double dgfr; - double dalq; - double dflo; -}; - -VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation& rhs); -VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation& rhs); -VFPEvaluation operator*(double lhs, const VFPEvaluation& rhs); - -/** - * Helper function which interpolates data using the indices etc. given in the inputs. - */ -VFPEvaluation interpolate(const VFPProdTable& table, - const InterpData& flo_i, - const InterpData& thp_i, - const InterpData& wfr_i, - const InterpData& gfr_i, - const InterpData& alq_i); - -/** - * This basically models interpolate(VFPProdTable::array_type, ...) - * which performs 5D interpolation, but here for the 2D case only - */ -VFPEvaluation interpolate(const VFPInjTable& table, - const InterpData& flo_i, - const InterpData& thp_i); - -VFPEvaluation bhp(const VFPProdTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp, - const double alq, - const double explicit_wfr, - const double explicit_gfr, - const bool use_vfpexplicit); - -VFPEvaluation bhp(const VFPInjTable& table, - const double aqua, - const double liquid, - const double vapour, - const double thp); - - /** * Returns the table from the map if found, or throws an exception */ @@ -164,54 +127,95 @@ bool hasTable(const std::map>& tables, int return (entry != tables.end() ); } - /** * Returns the type variable for FLO/GFR/WFR for production tables */ template TYPE getType(const TABLE& table); +} -/** - * This function finds the value of THP given a specific BHP. - * Essentially: - * Given the function f(thp_array(x)) = bhp_array(x), which is piecewise linear, - * find thp so that f(thp) = bhp. - */ -double findTHP(const std::vector& bhp_array, - const std::vector& thp_array, - double bhp); +template +class VFPHelpers { +public: + /** + * Helper function to find indices etc. for linear interpolation and extrapolation + * @param value_in Value to find in values + * @param values Sorted list of values to search for value in. + * @return Data required to find the interpolated value + */ + static detail::InterpData findInterpData(const Scalar value_in, + const std::vector& values); -/** -* Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq -*/ -std::pair -getMinimumBHPCoordinate(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq); + /** + * Helper function which interpolates data using the indices etc. given in the inputs. + */ + static detail::VFPEvaluation interpolate(const VFPProdTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i, + const detail::InterpData& wfr_i, + const detail::InterpData& gfr_i, + const detail::InterpData& alq_i); -/** -* Get (flo, bhp) at largest occuring stable vfp/ipr-intersection -* if it exists -*/ -std::optional> -intersectWithIPR(const VFPProdTable& table, - const double thp, - const double wfr, - const double gfr, - const double alq, - const double ipr_a, - const double ipr_b, - const std::function& adjust_bhp); + /** + * This basically models interpolate(VFPProdTable::array_type, ...) + * which performs 5D interpolation, but here for the 2D case only + */ + static detail::VFPEvaluation interpolate(const VFPInjTable& table, + const detail::InterpData& flo_i, + const detail::InterpData& thp_i); -} // namespace detail + static detail::VFPEvaluation bhp(const VFPProdTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp, + const Scalar alq, + const Scalar explicit_wfr, + const Scalar explicit_gfr, + const bool use_vfpexplicit); + static detail::VFPEvaluation bhp(const VFPInjTable& table, + const Scalar aqua, + const Scalar liquid, + const Scalar vapour, + const Scalar thp); + + /** + * This function finds the value of THP given a specific BHP. + * Essentially: + * Given the function f(thp_array(x)) = bhp_array(x), which is piecewise linear, + * find thp so that f(thp) = bhp. + */ + static Scalar findTHP(const std::vector& bhp_array, + const std::vector& thp_array, + Scalar bhp); + + /** + * Get (flo, bhp) at minimum bhp for given thp,wfr,gfr,alq + */ + static std::pair + getMinimumBHPCoordinate(const VFPProdTable& table, + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq); + + /** + * Get (flo, bhp) at largest occuring stable vfp/ipr-intersection + * if it exists + */ + static std::optional> + intersectWithIPR(const VFPProdTable& table, + const Scalar thp, + const Scalar wfr, + const Scalar gfr, + const Scalar alq, + const Scalar ipr_a, + const Scalar ipr_b, + const std::function& adjust_bhp); +}; } // namespace - - - #endif /* OPM_AUTODIFF_VFPHELPERS_HPP_ */ diff --git a/opm/simulators/wells/VFPInjProperties.cpp b/opm/simulators/wells/VFPInjProperties.cpp index 2da1e3749..0ba7be123 100644 --- a/opm/simulators/wells/VFPInjProperties.cpp +++ b/opm/simulators/wells/VFPInjProperties.cpp @@ -41,7 +41,7 @@ bhp(const int table_id, { const VFPInjTable& table = detail::getTable(m_tables, table_id); - detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg); + detail::VFPEvaluation retval = VFPHelpers::bhp(table, aqua, liquid, vapour, thp_arg); return retval.value; } @@ -61,7 +61,7 @@ thp(const int table_id, return 0.; } - const std::vector thp_array = table.getTHPAxis(); + const auto thp_array = table.getTHPAxis(); int nthp = thp_array.size(); /** @@ -69,14 +69,14 @@ thp(const 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 = detail::findInterpData(flo, table.getFloAxis()); + const auto flo_i = VFPHelpers::findInterpData(flo, table.getFloAxis()); std::vector bhp_array(nthp); for (int i = 0; i < nthp; ++i) { - auto thp_i = detail::findInterpData(thp_array[i], thp_array); - bhp_array[i] = detail::interpolate(table, flo_i, thp_i).value; + auto thp_i = VFPHelpers::findInterpData(thp_array[i], thp_array); + bhp_array[i] = VFPHelpers::interpolate(table, flo_i, thp_i).value; } - return detail::findTHP(bhp_array, thp_array, bhp_arg); + return VFPHelpers::findTHP(bhp_array, thp_array, bhp_arg); } template @@ -115,10 +115,10 @@ EvalWell VFPInjProperties::bhp(const int table_id, //First, find the values to interpolate between //Value of FLO is negative in OPM for producers, but positive in VFP table - auto flo_i = detail::findInterpData(flo.value(), table.getFloAxis()); - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant + const auto flo_i = VFPHelpers::findInterpData(flo.value(), table.getFloAxis()); + const auto thp_i = VFPHelpers::findInterpData(thp, table.getTHPAxis()); // assume constant - detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i); + detail::VFPEvaluation bhp_val = VFPHelpers::interpolate(table, flo_i, thp_i); bhp = bhp_val.dflo * flo; bhp.setValue(bhp_val.value); // thp is assumed constant i.e. diff --git a/opm/simulators/wells/VFPProdProperties.cpp b/opm/simulators/wells/VFPProdProperties.cpp index 8114d5b03..1ddec2d7c 100644 --- a/opm/simulators/wells/VFPProdProperties.cpp +++ b/opm/simulators/wells/VFPProdProperties.cpp @@ -67,17 +67,17 @@ thp(const 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 = 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 = VFPHelpers::findInterpData( flo, table.getFloAxis()); + auto wfr_i = VFPHelpers::findInterpData( wfr, table.getWFRAxis()); + auto gfr_i = VFPHelpers::findInterpData( gfr, table.getGFRAxis()); + auto alq_i = VFPHelpers::findInterpData( alq, table.getALQAxis()); std::vector bhp_array(nthp); for (int i = 0; i < nthp; ++i) { - auto thp_i = detail::findInterpData(thp_array[i], thp_array); - bhp_array[i] = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; + auto thp_i = VFPHelpers::findInterpData(thp_array[i], thp_array); + bhp_array[i] = VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i).value; } - return detail::findTHP(bhp_array, thp_array, bhp_arg); + return VFPHelpers::findTHP(bhp_array, thp_array, bhp_arg); } template @@ -94,7 +94,9 @@ bhp(const int table_id, { const VFPProdTable& table = detail::getTable(m_tables, table_id); - detail::VFPEvaluation retval = detail::bhp(table, aqua, liquid, vapour, thp_arg, alq, explicit_wfr,explicit_gfr, use_expvfp); + detail::VFPEvaluation retval = VFPHelpers::bhp(table, aqua, liquid, vapour, + thp_arg, alq, explicit_wfr, + explicit_gfr, use_expvfp); return retval.value; } @@ -124,16 +126,16 @@ bhpwithflo(const std::vector& flos, { // Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); - const auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant - const auto wfr_i = detail::findInterpData( wfr, table.getWFRAxis()); - const auto gfr_i = detail::findInterpData( gfr, table.getGFRAxis()); - const auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant + const auto thp_i = VFPHelpers::findInterpData( thp, table.getTHPAxis()); // assume constant + const auto wfr_i = VFPHelpers::findInterpData( wfr, table.getWFRAxis()); + const auto gfr_i = VFPHelpers::findInterpData( gfr, table.getGFRAxis()); + const auto alq_i = VFPHelpers::findInterpData( alq, table.getALQAxis()); //assume constant std::vector bhps(flos.size(), 0.); for (std::size_t i = 0; i < flos.size(); ++i) { // Value of FLO is negative in OPM for producers, but positive in VFP table - const auto flo_i = detail::findInterpData(-flos[i], table.getFloAxis()); - const detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + const auto flo_i = VFPHelpers::findInterpData(-flos[i], table.getFloAxis()); + const detail::VFPEvaluation bhp_val = VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); // TODO: this kind of breaks the conventions for the functions here by putting dp within the function bhps[i] = bhp_val.value - dp; @@ -152,7 +154,7 @@ minimumBHP(const int table_id, { // Get the table const VFPProdTable& table = detail::getTable(m_tables, table_id); - const auto retval = detail::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq); + const auto retval = VFPHelpers::getMinimumBHPCoordinate(table, thp, wfr, gfr, alq); // returned pair is (flo, bhp) return retval.second; } @@ -191,13 +193,14 @@ bhp(const int table_id, //First, find the values to interpolate between //Value of FLO is negative in OPM for producers, but positive in VFP table - auto flo_i = detail::findInterpData(-flo.value(), table.getFloAxis()); - auto thp_i = detail::findInterpData( thp, table.getTHPAxis()); // assume constant - auto wfr_i = detail::findInterpData( wfr.value(), table.getWFRAxis()); - auto gfr_i = detail::findInterpData( gfr.value(), table.getGFRAxis()); - auto alq_i = detail::findInterpData( alq, table.getALQAxis()); //assume constant + auto flo_i = VFPHelpers::findInterpData(-flo.value(), table.getFloAxis()); + auto thp_i = VFPHelpers::findInterpData( thp, table.getTHPAxis()); // assume constant + auto wfr_i = VFPHelpers::findInterpData( wfr.value(), table.getWFRAxis()); + auto gfr_i = VFPHelpers::findInterpData( gfr.value(), table.getGFRAxis()); + auto alq_i = VFPHelpers::findInterpData( alq, table.getALQAxis()); //assume constant - detail::VFPEvaluation bhp_val = detail::interpolate(table, flo_i, thp_i, wfr_i, gfr_i, alq_i); + detail::VFPEvaluation bhp_val = VFPHelpers::interpolate(table, flo_i, thp_i, wfr_i, + gfr_i, alq_i); bhp = (bhp_val.dwfr * wfr) + (bhp_val.dgfr * gfr) - (std::max(0.0, bhp_val.dflo) * flo); bhp.setValue(bhp_val.value); diff --git a/opm/simulators/wells/WellBhpThpCalculator.cpp b/opm/simulators/wells/WellBhpThpCalculator.cpp index 14a4cabaf..b90ddf8e9 100644 --- a/opm/simulators/wells/WellBhpThpCalculator.cpp +++ b/opm/simulators/wells/WellBhpThpCalculator.cpp @@ -919,7 +919,9 @@ isStableSolution(const WellState& well_state, const auto& table = well_.vfpProperties()->getProd()->getTable(controls.vfp_table_number); const bool use_vfpexplicit = well_.useVfpExplicit(); - detail::VFPEvaluation bhp = detail::bhp(table, aqua, liquid, vapour, thp, well_.getALQ(well_state), wfr, gfr, use_vfpexplicit); + auto bhp = VFPHelpers::bhp(table, aqua, liquid, vapour, thp, + well_.getALQ(well_state), wfr, gfr, + use_vfpexplicit); if (bhp.dflo >= 0) { return true; @@ -964,7 +966,10 @@ estimateStableBhp(const WellState& well_state, auto bhp_adjusted = [this, &thp, &dp_hydro](const Scalar bhp) { return bhp - dp_hydro + getVfpBhpAdjustment(bhp, thp); }; - const auto retval = detail::intersectWithIPR(table, thp, wfr, gfr, well_.getALQ(well_state), ipr.first, ipr.second, bhp_adjusted); + const auto retval = VFPHelpers::intersectWithIPR(table, thp, wfr, gfr, + well_.getALQ(well_state), + ipr.first, ipr.second, + bhp_adjusted); if (retval.has_value()) { // returned pair is (flo, bhp) return retval.value().second; diff --git a/tests/test_vfpproperties.cpp b/tests/test_vfpproperties.cpp index 44e17b047..fac8740a7 100644 --- a/tests/test_vfpproperties.cpp +++ b/tests/test_vfpproperties.cpp @@ -68,12 +68,12 @@ BOOST_AUTO_TEST_CASE(findInterpData) double first = 1; double last = 15; - Opm::detail::InterpData eval0 = Opm::detail::findInterpData(exact, values); - 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); - Opm::detail::InterpData eval4 = Opm::detail::findInterpData(first, values); - Opm::detail::InterpData eval5 = Opm::detail::findInterpData(last, values); + auto eval0 = Opm::VFPHelpers::findInterpData(exact, values); + auto eval1 = Opm::VFPHelpers::findInterpData(interpolate, values); + auto eval2 = Opm::VFPHelpers::findInterpData(extrapolate_left, values); + auto eval3 = Opm::VFPHelpers::findInterpData(extrapolate_right, values); + auto eval4 = Opm::VFPHelpers::findInterpData(first, values); + auto eval5 = Opm::VFPHelpers::findInterpData(last, values); BOOST_CHECK_EQUAL(eval0.ind_[0], 2); BOOST_CHECK_EQUAL(eval0.ind_[1], 3); @@ -109,7 +109,7 @@ BOOST_AUTO_TEST_SUITE_END() // HelperTests * values data is given at */ struct TrivialFixture { - typedef Opm::detail::VFPEvaluation VFPEvaluation; + typedef Opm::detail::VFPEvaluation VFPEvaluation; TrivialFixture() : table_ids(1, 1), thp_axis{0.0, 1.0},