From 9bfa39224fd23806c0774863f371962360c0d25f Mon Sep 17 00:00:00 2001 From: Kai Bao Date: Wed, 14 Nov 2018 15:48:38 +0100 Subject: [PATCH] re-organize the interface of function findIntersectionForBhp() hopefully, make it easier to use. and also, there is no mistake introduced. --- opm/autodiff/VFPHelpers.hpp | 58 +++++++++++++----------------- opm/autodiff/VFPProdProperties.cpp | 11 ++++-- 2 files changed, 33 insertions(+), 36 deletions(-) diff --git a/opm/autodiff/VFPHelpers.hpp b/opm/autodiff/VFPHelpers.hpp index 9c2dcfc3e..a2b7562e3 100644 --- a/opm/autodiff/VFPHelpers.hpp +++ b/opm/autodiff/VFPHelpers.hpp @@ -765,7 +765,7 @@ inline double findTHP( // a data type use to do the intersection calculation to get the intial bhp under THP control -struct DataPoint { +struct RateBhpPair { double rate; double bhp; }; @@ -773,7 +773,7 @@ struct DataPoint { // looking for a intersection point a line segment and a line, they are both defined with two points // it is copied from #include , which should be removed since it is only required by the lagacy polymer -inline bool findIntersection(const std::array& line_segment, const std::array& line, double& bhp) { +inline bool findIntersection(const std::array& line_segment, const std::array& line, double& bhp) { const double x1 = line_segment[0].rate; const double y1 = line_segment[0].bhp; const double x2 = line_segment[1].rate; @@ -802,41 +802,41 @@ inline bool findIntersection(const std::array& line_segment, const } // calculating the BHP from thp through the intersection of VFP curves and inflow performance relationship -inline bool findIntersectionForBhp(const std::vector&rate_samples, - const std::vector&bhp_samples, - const double flo_rate1, - const double flo_rate2, - const double bhp1, - const double bhp2, +inline bool findIntersectionForBhp(const std::vector& ratebhp_samples, + const std::array& ratebhp_twopoints_ipr, double& obtained_bhp) { - // there possibly two intersection point, then we choose the bigger one - // we choose the bigger one, then it will be the later one in the rate_samples + // there possibly two intersection point, then we choose the one corresponding with the bigger rate - const size_t num_samples = rate_samples.size(); - assert(num_samples == bhp_samples.size()); + const double bhp1 = ratebhp_twopoints_ipr[0].bhp; + const double rate1 = ratebhp_twopoints_ipr[0].rate; - assert(flo_rate1 != flo_rate2); + const double bhp2 = ratebhp_twopoints_ipr[1].bhp; + const double rate2 = ratebhp_twopoints_ipr[1].rate; + + assert(rate1 != rate2); + + const double line_slope = (bhp2 - bhp1) / (rate2 - rate1); - const double line_slope = (bhp2 - bhp1) / (flo_rate2 - flo_rate1); // line equation will be // bhp - bhp1 - line_slope * (flo_rate - flo_rate1) = 0 auto flambda = [&](const double flo_rate, const double bhp) { - return bhp - bhp1 - line_slope * (flo_rate - flo_rate1); + return bhp - bhp1 - line_slope * (flo_rate - rate1); }; int number_intersection_found = 0; int index_segment = 0; // the intersection segment that intersection happens - for (size_t i = 0; i < rate_samples.size() - 1; ++i) { - const double temp1 = flambda(rate_samples[i], bhp_samples[i]); - const double temp2 = flambda(rate_samples[i+1], bhp_samples[i+1]); + const size_t num_samples = ratebhp_samples.size(); + for (size_t i = 0; i < num_samples - 1; ++i) { + const double temp1 = flambda(ratebhp_samples[i].rate, ratebhp_samples[i].bhp); + const double temp2 = flambda(ratebhp_samples[i+1].rate, ratebhp_samples[i+1].bhp); if (temp1 * temp2 <= 0.) { // intersection happens // in theory there should be maximum two intersection points // while considering the situation == 0. here, we might find more - // we always use the last one, which is the one has the biggest rate + // we always use the last one, which is the one corresponds to the biggest rate, + // which we assume is the more stable one ++number_intersection_found; index_segment = i; - std::cout << " temp1 " << temp1 << " temp2 " << temp2 << std::endl; } } @@ -844,22 +844,12 @@ inline bool findIntersectionForBhp(const std::vector&rate_samples, return false; } - // then we need to calculate the intersection point - const std::array line_segment{ DataPoint{rate_samples[index_segment], bhp_samples[index_segment]}, - DataPoint{rate_samples[index_segment + 1], bhp_samples[index_segment + 1]} }; + // then we pick the segment from the VFP curve to do the line intersection calculation + const std::array line_segment{ ratebhp_samples[index_segment], ratebhp_samples[index_segment + 1] }; - const std::array line { DataPoint{flo_rate1, bhp1}, - DataPoint{flo_rate2, bhp2} }; + const bool intersection_found = findIntersection(line_segment, ratebhp_twopoints_ipr, obtained_bhp); - const bool inter_section_found = findIntersection(line_segment, line, obtained_bhp); - - if (inter_section_found) { - std::cout << " found a bhp is " << obtained_bhp << std::endl; - return true; - } else { - std::cout << " did not find the intersection point " << std::endl; - return false; - } + return intersection_found; } diff --git a/opm/autodiff/VFPProdProperties.cpp b/opm/autodiff/VFPProdProperties.cpp index 50b1847c6..59bd406c3 100644 --- a/opm/autodiff/VFPProdProperties.cpp +++ b/opm/autodiff/VFPProdProperties.cpp @@ -220,9 +220,16 @@ calculateBhpWithTHPTarget(const std::vector& ipr_a, // get the bhp sampling values based on the flo sample values const std::vector bhp_flo_samples = bhpwithflo(flo_samples, thp_table_id, wfr, gfr, thp_limit, alq, dp); + std::vector ratebhp_samples; + for (size_t i = 0; i < flo_samples.size(); ++i) { + ratebhp_samples.push_back( detail::RateBhpPair{flo_samples[i], bhp_flo_samples[i]} ); + } + + const std::array ratebhp_twopoints_ipr {detail::RateBhpPair{flo_bhp_middle, bhp_middle}, + detail::RateBhpPair{flo_bhp_limit, bhp_limit} }; + double obtain_bhp = 0.; - const bool obtain_solution_with_thp_limit = detail::findIntersectionForBhp(flo_samples, bhp_flo_samples, - flo_bhp_middle, flo_bhp_limit, bhp_middle, bhp_limit, obtain_bhp); + const bool obtain_solution_with_thp_limit = detail::findIntersectionForBhp(ratebhp_samples, ratebhp_twopoints_ipr, obtain_bhp); // \Note: assuming not that negative BHP does not make sense if (obtain_solution_with_thp_limit && obtain_bhp > 0.) {