00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
00022 #define OPM_AUTODIFF_VFPHELPERS_HPP_
00023
00024
00025 #include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
00026 #include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
00027 #include <opm/autodiff/AutoDiffHelpers.hpp>
00028 #include <opm/material/densead/Math.hpp>
00029 #include <opm/material/densead/Evaluation.hpp>
00030
00034 namespace Opm {
00035 namespace detail {
00036
00037
00038 typedef AutoDiffBlock<double> ADB;
00039
00040
00044 inline double zeroIfNan(const double& value) {
00045 return (std::isnan(value)) ? 0.0 : value;
00046 }
00047
00051 template <class EvalWell>
00052 inline double zeroIfNan(const EvalWell& value) {
00053 return (std::isnan(value.value())) ? 0.0 : value.value();
00054 }
00055
00056
00057
00061 inline ADB zeroIfNan(const ADB& values) {
00062 Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
00063
00064 const ADB::V z = ADB::V::Zero(values.size());
00065 const ADB zero = ADB::constant(z, values.blockPattern());
00066
00067 ADB retval = not_nan_selector.select(values, zero);
00068 return retval;
00069 }
00070
00071
00072
00073
00074
00080 template <typename T>
00081 static T getFlo(const T& aqua, const T& liquid, const T& vapour,
00082 const VFPProdTable::FLO_TYPE& type) {
00083 switch (type) {
00084 case VFPProdTable::FLO_OIL:
00085
00086 return liquid;
00087 case VFPProdTable::FLO_LIQ:
00088
00089 return aqua + liquid;
00090 case VFPProdTable::FLO_GAS:
00091
00092 return vapour;
00093 case VFPProdTable::FLO_INVALID:
00094 default:
00095 OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
00096 }
00097 }
00098
00099
00100
00106 template <typename T>
00107 static T getFlo(const T& aqua, const T& liquid, const T& vapour,
00108 const VFPInjTable::FLO_TYPE& type) {
00109 switch (type) {
00110 case VFPInjTable::FLO_OIL:
00111
00112 return liquid;
00113 case VFPInjTable::FLO_WAT:
00114
00115 return aqua;
00116 case VFPInjTable::FLO_GAS:
00117
00118 return vapour;
00119 case VFPInjTable::FLO_INVALID:
00120 default:
00121 OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
00122 }
00123 }
00124
00125
00126
00127
00128
00129
00130
00135 template <typename T>
00136 static T getWFR(const T& aqua, const T& liquid, const T& vapour,
00137 const VFPProdTable::WFR_TYPE& type) {
00138 switch(type) {
00139 case VFPProdTable::WFR_WOR: {
00140
00141 T wor = aqua / liquid;
00142 return zeroIfNan(wor);
00143 }
00144 case VFPProdTable::WFR_WCT:
00145
00146 return zeroIfNan(aqua / (aqua + liquid));
00147 case VFPProdTable::WFR_WGR:
00148
00149 return zeroIfNan(aqua / vapour);
00150 case VFPProdTable::WFR_INVALID:
00151 default:
00152 OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'");
00153 }
00154 }
00155
00156
00157
00158
00159
00160
00165 template <typename T>
00166 static T getGFR(const T& aqua, const T& liquid, const T& vapour,
00167 const VFPProdTable::GFR_TYPE& type) {
00168 switch(type) {
00169 case VFPProdTable::GFR_GOR:
00170
00171 return zeroIfNan(vapour / liquid);
00172 case VFPProdTable::GFR_GLR:
00173
00174 return zeroIfNan(vapour / (liquid + aqua));
00175 case VFPProdTable::GFR_OGR:
00176
00177 return zeroIfNan(liquid / vapour);
00178 case VFPProdTable::GFR_INVALID:
00179 default:
00180 OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
00181 }
00182 }
00183
00184
00185
00186
00187
00188
00192 struct InterpData {
00193 InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {}
00194 int ind_[2];
00195 double inv_dist_;
00196 double factor_;
00197 };
00198
00199
00200
00201
00202
00203
00210 inline InterpData findInterpData(const double& value, const std::vector<double>& values) {
00211 InterpData retval;
00212
00213 const int nvalues = values.size();
00214
00215
00216 if (values.size() == 1) {
00217 retval.ind_[0] = 0;
00218 retval.ind_[1] = 0;
00219 retval.inv_dist_ = 0.0;
00220 retval.factor_ = 0.0;
00221 }
00222
00223 else {
00224
00225 if (value < values.front()) {
00226 retval.ind_[0] = 0;
00227 retval.ind_[1] = 1;
00228 }
00229
00230 else if (value >= values.back()) {
00231 retval.ind_[0] = nvalues-2;
00232 retval.ind_[1] = nvalues-1;
00233 }
00234 else {
00235
00236 for (int i=1; i<nvalues; ++i) {
00237 if (values[i] >= value) {
00238 retval.ind_[0] = i-1;
00239 retval.ind_[1] = i;
00240 break;
00241 }
00242 }
00243 }
00244
00245 const double start = values[retval.ind_[0]];
00246 const double end = values[retval.ind_[1]];
00247
00248
00249 if (end > start) {
00250
00251
00252 retval.inv_dist_ = 1.0 / (end-start);
00253 retval.factor_ = (value-start) * retval.inv_dist_;
00254 }
00255 else {
00256 retval.inv_dist_ = 0.0;
00257 retval.factor_ = 0.0;
00258 }
00259 }
00260
00261 return retval;
00262 }
00263
00264
00265
00266
00267
00268
00272 struct VFPEvaluation {
00273 VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
00274 double value;
00275 double dthp;
00276 double dwfr;
00277 double dgfr;
00278 double dalq;
00279 double dflo;
00280 };
00281
00282 inline VFPEvaluation operator+(
00283 VFPEvaluation lhs,
00284 const VFPEvaluation& rhs) {
00285 lhs.value += rhs.value;
00286 lhs.dthp += rhs.dthp;
00287 lhs.dwfr += rhs.dwfr;
00288 lhs.dgfr += rhs.dgfr;
00289 lhs.dalq += rhs.dalq;
00290 lhs.dflo += rhs.dflo;
00291 return lhs;
00292 }
00293
00294 inline VFPEvaluation operator-(
00295 VFPEvaluation lhs,
00296 const VFPEvaluation& rhs) {
00297 lhs.value -= rhs.value;
00298 lhs.dthp -= rhs.dthp;
00299 lhs.dwfr -= rhs.dwfr;
00300 lhs.dgfr -= rhs.dgfr;
00301 lhs.dalq -= rhs.dalq;
00302 lhs.dflo -= rhs.dflo;
00303 return lhs;
00304 }
00305
00306 inline VFPEvaluation operator*(
00307 double lhs,
00308 const VFPEvaluation& rhs) {
00309 VFPEvaluation retval;
00310 retval.value = rhs.value * lhs;
00311 retval.dthp = rhs.dthp * lhs;
00312 retval.dwfr = rhs.dwfr * lhs;
00313 retval.dgfr = rhs.dgfr * lhs;
00314 retval.dalq = rhs.dalq * lhs;
00315 retval.dflo = rhs.dflo * lhs;
00316 return retval;
00317 }
00318
00319
00320
00321
00322
00323
00327 #ifdef __GNUC__
00328 #pragma GCC push_options
00329 #pragma GCC optimize ("unroll-loops")
00330 #endif
00331
00332
00333 inline VFPEvaluation interpolate(
00334 const VFPProdTable::array_type& array,
00335 const InterpData& flo_i,
00336 const InterpData& thp_i,
00337 const InterpData& wfr_i,
00338 const InterpData& gfr_i,
00339 const InterpData& alq_i) {
00340
00341
00342 VFPEvaluation nn[2][2][2][2][2];
00343
00344
00345
00346
00347
00348
00349
00350 for (int t=0; t<=1; ++t) {
00351 for (int w=0; w<=1; ++w) {
00352 for (int g=0; g<=1; ++g) {
00353 for (int a=0; a<=1; ++a) {
00354 for (int f=0; f<=1; ++f) {
00355
00356 const int ti = thp_i.ind_[t];
00357 const int wi = wfr_i.ind_[w];
00358 const int gi = gfr_i.ind_[g];
00359 const int ai = alq_i.ind_[a];
00360 const int fi = flo_i.ind_[f];
00361
00362
00363 nn[t][w][g][a][f].value = array[ti][wi][gi][ai][fi];
00364 }
00365 }
00366 }
00367 }
00368 }
00369
00370
00371
00372
00373 for (int i=0; i<=1; ++i) {
00374 for (int j=0; j<=1; ++j) {
00375 for (int k=0; k<=1; ++k) {
00376 for (int l=0; l<=1; ++l) {
00377 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_;
00378 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_;
00379 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_;
00380 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_;
00381 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_;
00382
00383 nn[1][i][j][k][l].dthp = nn[0][i][j][k][l].dthp;
00384 nn[i][1][j][k][l].dwfr = nn[i][0][j][k][l].dwfr;
00385 nn[i][j][1][k][l].dgfr = nn[i][j][0][k][l].dgfr;
00386 nn[i][j][k][1][l].dalq = nn[i][j][k][0][l].dalq;
00387 nn[i][j][k][l][1].dflo = nn[i][j][k][l][0].dflo;
00388 }
00389 }
00390 }
00391 }
00392
00393 double t1, t2;
00394
00395
00396
00397
00398
00399 t2 = flo_i.factor_;
00400 t1 = (1.0-t2);
00401 for (int t=0; t<=1; ++t) {
00402 for (int w=0; w<=1; ++w) {
00403 for (int g=0; g<=1; ++g) {
00404 for (int a=0; a<=1; ++a) {
00405 nn[t][w][g][a][0] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1];
00406 }
00407 }
00408 }
00409 }
00410
00411 t2 = alq_i.factor_;
00412 t1 = (1.0-t2);
00413 for (int t=0; t<=1; ++t) {
00414 for (int w=0; w<=1; ++w) {
00415 for (int g=0; g<=1; ++g) {
00416 nn[t][w][g][0][0] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0];
00417 }
00418 }
00419 }
00420
00421 t2 = gfr_i.factor_;
00422 t1 = (1.0-t2);
00423 for (int t=0; t<=1; ++t) {
00424 for (int w=0; w<=1; ++w) {
00425 nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0];
00426 }
00427 }
00428
00429 t2 = wfr_i.factor_;
00430 t1 = (1.0-t2);
00431 for (int t=0; t<=1; ++t) {
00432 nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0];
00433 }
00434
00435 t2 = thp_i.factor_;
00436 t1 = (1.0-t2);
00437 nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0];
00438
00439 return nn[0][0][0][0][0];
00440 }
00441
00442
00443
00444
00445
00450 inline VFPEvaluation interpolate(
00451 const VFPInjTable::array_type& array,
00452 const InterpData& flo_i,
00453 const InterpData& thp_i) {
00454
00455
00456 VFPEvaluation nn[2][2];
00457
00458
00459
00460
00461 for (int t=0; t<=1; ++t) {
00462 for (int f=0; f<=1; ++f) {
00463
00464 const int ti = thp_i.ind_[t];
00465 const int fi = flo_i.ind_[f];
00466
00467
00468 nn[t][f].value = array[ti][fi];
00469 }
00470 }
00471
00472
00473
00474
00475 for (int i=0; i<=1; ++i) {
00476 nn[0][i].dthp = (nn[1][i].value - nn[0][i].value) * thp_i.inv_dist_;
00477 nn[i][0].dwfr = -1e100;
00478 nn[i][0].dgfr = -1e100;
00479 nn[i][0].dalq = -1e100;
00480 nn[i][0].dflo = (nn[i][1].value - nn[i][0].value) * flo_i.inv_dist_;
00481
00482 nn[1][i].dthp = nn[0][i].dthp;
00483 nn[i][1].dwfr = nn[i][0].dwfr;
00484 nn[i][1].dgfr = nn[i][0].dgfr;
00485 nn[i][1].dalq = nn[i][0].dalq;
00486 nn[i][1].dflo = nn[i][0].dflo;
00487 }
00488
00489 double t1, t2;
00490
00491
00492 t2 = flo_i.factor_;
00493 t1 = (1.0-t2);
00494 nn[0][0] = t1*nn[0][0] + t2*nn[0][1];
00495 nn[1][0] = t1*nn[1][0] + t2*nn[1][1];
00496
00497
00498 t2 = thp_i.factor_;
00499 t1 = (1.0-t2);
00500 nn[0][0] = t1*nn[0][0] + t2*nn[1][0];
00501
00502 return nn[0][0];
00503 }
00504
00505
00506
00507
00508 #ifdef __GNUC__
00509 #pragma GCC pop_options //unroll loops
00510 #endif
00511
00512
00513
00514
00515
00516 inline VFPEvaluation bhp(const VFPProdTable* table,
00517 const double& aqua,
00518 const double& liquid,
00519 const double& vapour,
00520 const double& thp,
00521 const double& alq) {
00522
00523 double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
00524 double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
00525 double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
00526
00527
00528
00529 auto flo_i = detail::findInterpData(-flo, table->getFloAxis());
00530 auto thp_i = detail::findInterpData( thp, table->getTHPAxis());
00531 auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis());
00532 auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis());
00533 auto alq_i = detail::findInterpData( alq, table->getALQAxis());
00534
00535 detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
00536
00537 return retval;
00538 }
00539
00540
00541
00542
00543
00544 inline VFPEvaluation bhp(const VFPInjTable* table,
00545 const double& aqua,
00546 const double& liquid,
00547 const double& vapour,
00548 const double& thp) {
00549
00550 double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
00551
00552
00553 auto flo_i = detail::findInterpData(flo, table->getFloAxis());
00554 auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
00555
00556
00557 detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i);
00558
00559 return retval;
00560 }
00561
00562
00563
00564
00565
00566
00567
00568
00572 template <typename T>
00573 const T* getTable(const std::map<int, T*> tables, int table_id) {
00574 auto entry = tables.find(table_id);
00575 if (entry == tables.end()) {
00576 OPM_THROW(std::invalid_argument, "Nonexistent table " << table_id << " referenced.");
00577 }
00578 else {
00579 return entry->second;
00580 }
00581 }
00582
00583
00584
00585
00586
00587
00588
00589
00590
00591
00595 inline void extendBlockPattern(const ADB& x, std::vector<int>& block_pattern) {
00596 std::vector<int> x_block_pattern = x.blockPattern();
00597
00598 if (x_block_pattern.empty()) {
00599 return;
00600 }
00601 else {
00602 if (block_pattern.empty()) {
00603 block_pattern = x_block_pattern;
00604 return;
00605 }
00606 else {
00607 if (x_block_pattern != block_pattern) {
00608 OPM_THROW(std::logic_error, "Block patterns do not match");
00609 }
00610 }
00611 }
00612 }
00613
00617 inline std::vector<int> commonBlockPattern(
00618 const ADB& x1,
00619 const ADB& x2,
00620 const ADB& x3,
00621 const ADB& x4) {
00622 std::vector<int> block_pattern;
00623
00624 extendBlockPattern(x1, block_pattern);
00625 extendBlockPattern(x2, block_pattern);
00626 extendBlockPattern(x3, block_pattern);
00627 extendBlockPattern(x4, block_pattern);
00628
00629 return block_pattern;
00630 }
00631
00632 inline std::vector<int> commonBlockPattern(
00633 const ADB& x1,
00634 const ADB& x2,
00635 const ADB& x3,
00636 const ADB& x4,
00637 const ADB& x5) {
00638 std::vector<int> block_pattern = commonBlockPattern(x1, x2, x3, x4);
00639 extendBlockPattern(x5, block_pattern);
00640
00641 return block_pattern;
00642 }
00643
00644
00645
00646
00647
00648
00649
00650
00651
00652
00653
00654
00658 template <typename TYPE, typename TABLE>
00659 TYPE getType(const TABLE* table);
00660
00661 template <>
00662 inline
00663 VFPProdTable::FLO_TYPE getType(const VFPProdTable* table) {
00664 return table->getFloType();
00665 }
00666
00667 template <>
00668 inline
00669 VFPProdTable::WFR_TYPE getType(const VFPProdTable* table) {
00670 return table->getWFRType();
00671 }
00672
00673 template <>
00674 inline
00675 VFPProdTable::GFR_TYPE getType(const VFPProdTable* table) {
00676 return table->getGFRType();
00677 }
00678
00679
00683 template <>
00684 inline
00685 VFPInjTable::FLO_TYPE getType(const VFPInjTable* table) {
00686 return table->getFloType();
00687 }
00688
00689
00690
00691
00695 template <typename TYPE>
00696 ADB getValue(
00697 const ADB& aqua,
00698 const ADB& liquid,
00699 const ADB& vapour, TYPE type);
00700
00701 template <>
00702 inline
00703 ADB getValue(
00704 const ADB& aqua,
00705 const ADB& liquid,
00706 const ADB& vapour,
00707 VFPProdTable::FLO_TYPE type) {
00708 return detail::getFlo(aqua, liquid, vapour, type);
00709 }
00710
00711 template <>
00712 inline
00713 ADB getValue(
00714 const ADB& aqua,
00715 const ADB& liquid,
00716 const ADB& vapour,
00717 VFPProdTable::WFR_TYPE type) {
00718 return detail::getWFR(aqua, liquid, vapour, type);
00719 }
00720
00721 template <>
00722 inline
00723 ADB getValue(
00724 const ADB& aqua,
00725 const ADB& liquid,
00726 const ADB& vapour,
00727 VFPProdTable::GFR_TYPE type) {
00728 return detail::getGFR(aqua, liquid, vapour, type);
00729 }
00730
00731 template <>
00732 inline
00733 ADB getValue(
00734 const ADB& aqua,
00735 const ADB& liquid,
00736 const ADB& vapour,
00737 VFPInjTable::FLO_TYPE type) {
00738 return detail::getFlo(aqua, liquid, vapour, type);
00739 }
00740
00748 template <typename TYPE, typename TABLE>
00749 ADB combineADBVars(const std::vector<const TABLE*>& well_tables,
00750 const ADB& aqua,
00751 const ADB& liquid,
00752 const ADB& vapour) {
00753
00754 const int num_wells = static_cast<int>(well_tables.size());
00755 assert(aqua.size() == num_wells);
00756 assert(liquid.size() == num_wells);
00757 assert(vapour.size() == num_wells);
00758
00759
00760 std::map<TYPE, ADB> map;
00761
00762
00763 std::map<TYPE, std::vector<int> > elems;
00764
00765
00766
00767 for (int i=0; i<num_wells; ++i) {
00768 const TABLE* table = well_tables[i];
00769
00770
00771 if (table != NULL) {
00772 TYPE type = getType<TYPE>(table);
00773
00774
00775
00776 if (map.find(type) == map.end()) {
00777 map.insert(std::pair<TYPE, ADB>(
00778 type,
00779 detail::getValue<TYPE>(aqua, liquid, vapour, type)
00780 ));
00781 }
00782
00783
00784 elems[type].push_back(i);
00785 }
00786 }
00787
00788
00789
00790 ADB retval = ADB::constant(ADB::V::Zero(num_wells));
00791 for (const auto& entry : elems) {
00792 const auto& key = entry.first;
00793 const auto& value = entry.second;
00794
00795
00796 assert(map.find(key) != map.end());
00797 const ADB& values = map.find(key)->second;
00798
00799
00800 const std::vector<int>& current = value;
00801
00802
00803 retval = retval + superset(subset(values, current), current, values.size());
00804 }
00805
00806 return retval;
00807 }
00808
00813 inline double findX(const double& x0,
00814 const double& x1,
00815 const double& y0,
00816 const double& y1,
00817 const double& y) {
00818 const double dx = x1 - x0;
00819 const double dy = y1 - y0;
00820
00828 double x = 0.0;
00829
00830 if (dy != 0.0) {
00831 x = x0 + (y-y0) * (dx/dy);
00832 }
00833 else {
00834 x = x1;
00835 }
00836
00837 return x;
00838 }
00839
00840
00841
00842
00843
00844
00845
00846
00847
00848
00855 inline double findTHP(
00856 const std::vector<double>& bhp_array,
00857 const std::vector<double>& thp_array,
00858 double bhp) {
00859 int nthp = thp_array.size();
00860
00861 double thp = -1e100;
00862
00863
00864 assert(std::is_sorted(thp_array.begin(), thp_array.end()));
00865
00872 if (std::is_sorted(bhp_array.begin(), bhp_array.end())) {
00873
00874 if (bhp <= bhp_array[0]) {
00875
00876 const double& x0 = thp_array[0];
00877 const double& x1 = thp_array[1];
00878 const double& y0 = bhp_array[0];
00879 const double& y1 = bhp_array[1];
00880 thp = detail::findX(x0, x1, y0, y1, bhp);
00881 }
00882
00883 else if (bhp > bhp_array[nthp-1]) {
00884
00885 const double& x0 = thp_array[nthp-2];
00886 const double& x1 = thp_array[nthp-1];
00887 const double& y0 = bhp_array[nthp-2];
00888 const double& y1 = bhp_array[nthp-1];
00889 thp = detail::findX(x0, x1, y0, y1, bhp);
00890 }
00891
00892 else {
00893
00894
00895
00896
00897
00898
00899 int i=0;
00900 bool found = false;
00901 for (; i<nthp-1; ++i) {
00902 const double& y0 = bhp_array[i ];
00903 const double& y1 = bhp_array[i+1];
00904
00905 if (y0 < bhp && bhp <= y1) {
00906 found = true;
00907 break;
00908 }
00909 }
00910
00911 assert(found == true);
00912 static_cast<void>(found);
00913
00914 const double& x0 = thp_array[i ];
00915 const double& x1 = thp_array[i+1];
00916 const double& y0 = bhp_array[i ];
00917 const double& y1 = bhp_array[i+1];
00918 thp = detail::findX(x0, x1, y0, y1, bhp);
00919 }
00920 }
00921
00922 else {
00923
00924
00925
00926 int i=0;
00927 bool found = false;
00928 for (; i<nthp-1; ++i) {
00929 const double& y0 = bhp_array[i ];
00930 const double& y1 = bhp_array[i+1];
00931
00932 if (y0 < bhp && bhp <= y1) {
00933 found = true;
00934 break;
00935 }
00936 }
00937 if (found) {
00938 const double& x0 = thp_array[i ];
00939 const double& x1 = thp_array[i+1];
00940 const double& y0 = bhp_array[i ];
00941 const double& y1 = bhp_array[i+1];
00942 thp = detail::findX(x0, x1, y0, y1, bhp);
00943 }
00944 else if (bhp <= bhp_array[0]) {
00945
00946 const double& x0 = thp_array[0];
00947 const double& x1 = thp_array[1];
00948 const double& y0 = bhp_array[0];
00949 const double& y1 = bhp_array[1];
00950 thp = detail::findX(x0, x1, y0, y1, bhp);
00951 }
00952
00953 else if (bhp > bhp_array[nthp-1]) {
00954
00955 const double& x0 = thp_array[nthp-2];
00956 const double& x1 = thp_array[nthp-1];
00957 const double& y0 = bhp_array[nthp-2];
00958 const double& y1 = bhp_array[nthp-1];
00959 thp = detail::findX(x0, x1, y0, y1, bhp);
00960 }
00961 else {
00962 OPM_THROW(std::logic_error, "Programmer error: Unable to find THP in THP array");
00963 }
00964 }
00965
00966 return thp;
00967 }
00968
00969
00970
00971
00972
00973
00974
00975 }
00976
00977
00978 }
00979
00980
00981
00982
00983 #endif