OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

re-organize the interface of function findIntersectionForBhp()

Browse Source 
hopefully, make it easier to use.

and also, there is no mistake introduced.
This commit is contained in:
Kai Bao 2018-11-14 15:48:38 +01:00
parent 8445c802c0
commit 9bfa39224f
2 changed files with 33 additions and 36 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

58
opm/autodiff/VFPHelpers.hpp
View File

@ -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 <opm/polymer/Point2D.hpp>, which should be removed since it is only required by the lagacy polymer
inline bool findIntersection(const std::array<DataPoint, 2>& line_segment, const std::array<DataPoint, 2>& line, double& bhp) {
inline bool findIntersection(const std::array<RateBhpPair, 2>& line_segment, const std::array<RateBhpPair, 2>& 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<DataPoint, 2>& line_segment, const
}
// calculating the BHP from thp through the intersection of VFP curves and inflow performance relationship
inline bool findIntersectionForBhp(const std::vector<double>&rate_samples,
const std::vector<double>&bhp_samples,
const double flo_rate1,
const double flo_rate2,
const double bhp1,
const double bhp2,
inline bool findIntersectionForBhp(const std::vector<RateBhpPair>& ratebhp_samples,
const std::array<RateBhpPair, 2>& 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<double>&rate_samples,
return false;
}
// then we need to calculate the intersection point
const std::array<DataPoint, 2> 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<RateBhpPair, 2> line_segment{ ratebhp_samples[index_segment], ratebhp_samples[index_segment + 1] };
const std::array<DataPoint, 2> 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;
}

11
opm/autodiff/VFPProdProperties.cpp
View File

@ -220,9 +220,16 @@ calculateBhpWithTHPTarget(const std::vector<double>& ipr_a,
// get the bhp sampling values based on the flo sample values
const std::vector<double> bhp_flo_samples = bhpwithflo(flo_samples, thp_table_id, wfr, gfr, thp_limit, alq, dp);
std::vector<detail::RateBhpPair> 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<detail::RateBhpPair, 2> 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.) {