21 #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
22 #define OPM_AUTODIFF_VFPHELPERS_HPP_
25 #include <opm/parser/eclipse/EclipseState/Tables/VFPProdTable.hpp>
26 #include <opm/parser/eclipse/EclipseState/Tables/VFPInjTable.hpp>
27 #include <opm/autodiff/AutoDiffHelpers.hpp>
28 #include <opm/material/densead/Math.hpp>
29 #include <opm/material/densead/Evaluation.hpp>
38 typedef AutoDiffBlock<double> ADB;
44 inline double zeroIfNan(
const double& value) {
45 return (std::isnan(value)) ? 0.0 : value;
51 template <
class EvalWell>
52 inline double zeroIfNan(
const EvalWell& value) {
53 return (std::isnan(value.value())) ? 0.0 : value.value();
61 inline ADB zeroIfNan(
const ADB& values) {
62 Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
64 const ADB::V z = ADB::V::Zero(values.size());
67 ADB retval = not_nan_selector.select(values, zero);
81 static T getFlo(
const T& aqua,
const T& liquid,
const T& vapour,
82 const VFPProdTable::FLO_TYPE& type) {
84 case VFPProdTable::FLO_OIL:
87 case VFPProdTable::FLO_LIQ:
90 case VFPProdTable::FLO_GAS:
93 case VFPProdTable::FLO_INVALID:
95 OPM_THROW(std::logic_error,
"Invalid FLO_TYPE: '" << type <<
"'");
106 template <
typename T>
107 static T getFlo(
const T& aqua,
const T& liquid,
const T& vapour,
108 const VFPInjTable::FLO_TYPE& type) {
110 case VFPInjTable::FLO_OIL:
113 case VFPInjTable::FLO_WAT:
116 case VFPInjTable::FLO_GAS:
119 case VFPInjTable::FLO_INVALID:
121 OPM_THROW(std::logic_error,
"Invalid FLO_TYPE: '" << type <<
"'");
135 template <
typename T>
136 static T getWFR(
const T& aqua,
const T& liquid,
const T& vapour,
137 const VFPProdTable::WFR_TYPE& type) {
139 case VFPProdTable::WFR_WOR: {
141 T wor = aqua / liquid;
142 return zeroIfNan(wor);
144 case VFPProdTable::WFR_WCT:
146 return zeroIfNan(aqua / (aqua + liquid));
147 case VFPProdTable::WFR_WGR:
149 return zeroIfNan(aqua / vapour);
150 case VFPProdTable::WFR_INVALID:
152 OPM_THROW(std::logic_error,
"Invalid WFR_TYPE: '" << type <<
"'");
165 template <
typename T>
166 static T getGFR(
const T& aqua,
const T& liquid,
const T& vapour,
167 const VFPProdTable::GFR_TYPE& type) {
169 case VFPProdTable::GFR_GOR:
171 return zeroIfNan(vapour / liquid);
172 case VFPProdTable::GFR_GLR:
174 return zeroIfNan(vapour / (liquid + aqua));
175 case VFPProdTable::GFR_OGR:
177 return zeroIfNan(liquid / vapour);
178 case VFPProdTable::GFR_INVALID:
180 OPM_THROW(std::logic_error,
"Invalid GFR_TYPE: '" << type <<
"'");
193 InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {}
210 inline InterpData findInterpData(
const double& value,
const std::vector<double>& values) {
213 const int nvalues = values.size();
216 if (values.size() == 1) {
219 retval.inv_dist_ = 0.0;
220 retval.factor_ = 0.0;
225 if (value < values.front()) {
230 else if (value >= values.back()) {
231 retval.ind_[0] = nvalues-2;
232 retval.ind_[1] = nvalues-1;
236 for (
int i=1; i<nvalues; ++i) {
237 if (values[i] >= value) {
238 retval.ind_[0] = i-1;
245 const double start = values[retval.ind_[0]];
246 const double end = values[retval.ind_[1]];
252 retval.inv_dist_ = 1.0 / (end-start);
253 retval.factor_ = (value-start) * retval.inv_dist_;
256 retval.inv_dist_ = 0.0;
257 retval.factor_ = 0.0;
273 VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
285 lhs.value += rhs.value;
286 lhs.dthp += rhs.dthp;
287 lhs.dwfr += rhs.dwfr;
288 lhs.dgfr += rhs.dgfr;
289 lhs.dalq += rhs.dalq;
290 lhs.dflo += rhs.dflo;
294 inline VFPEvaluation operator-(
296 const VFPEvaluation& rhs) {
297 lhs.value -= rhs.value;
298 lhs.dthp -= rhs.dthp;
299 lhs.dwfr -= rhs.dwfr;
300 lhs.dgfr -= rhs.dgfr;
301 lhs.dalq -= rhs.dalq;
302 lhs.dflo -= rhs.dflo;
306 inline VFPEvaluation operator*(
308 const VFPEvaluation& rhs) {
309 VFPEvaluation retval;
310 retval.value = rhs.value * lhs;
311 retval.dthp = rhs.dthp * lhs;
312 retval.dwfr = rhs.dwfr * lhs;
313 retval.dgfr = rhs.dgfr * lhs;
314 retval.dalq = rhs.dalq * lhs;
315 retval.dflo = rhs.dflo * lhs;
328 #pragma GCC push_options
329 #pragma GCC optimize ("unroll-loops")
333 inline VFPEvaluation interpolate(
334 const VFPProdTable::array_type& array,
335 const InterpData& flo_i,
336 const InterpData& thp_i,
337 const InterpData& wfr_i,
338 const InterpData& gfr_i,
339 const InterpData& alq_i) {
342 VFPEvaluation nn[2][2][2][2][2];
350 for (
int t=0; t<=1; ++t) {
351 for (
int w=0; w<=1; ++w) {
352 for (
int g=0; g<=1; ++g) {
353 for (
int a=0; a<=1; ++a) {
354 for (
int f=0; f<=1; ++f) {
356 const int ti = thp_i.ind_[t];
357 const int wi = wfr_i.ind_[w];
358 const int gi = gfr_i.ind_[g];
359 const int ai = alq_i.ind_[a];
360 const int fi = flo_i.ind_[f];
363 nn[t][w][g][a][f].value = array[ti][wi][gi][ai][fi];
373 for (
int i=0; i<=1; ++i) {
374 for (
int j=0; j<=1; ++j) {
375 for (
int k=0; k<=1; ++k) {
376 for (
int l=0; l<=1; ++l) {
377 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_;
378 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_;
379 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_;
380 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_;
381 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_;
383 nn[1][i][j][k][l].dthp = nn[0][i][j][k][l].dthp;
384 nn[i][1][j][k][l].dwfr = nn[i][0][j][k][l].dwfr;
385 nn[i][j][1][k][l].dgfr = nn[i][j][0][k][l].dgfr;
386 nn[i][j][k][1][l].dalq = nn[i][j][k][0][l].dalq;
387 nn[i][j][k][l][1].dflo = nn[i][j][k][l][0].dflo;
401 for (
int t=0; t<=1; ++t) {
402 for (
int w=0; w<=1; ++w) {
403 for (
int g=0; g<=1; ++g) {
404 for (
int a=0; a<=1; ++a) {
405 nn[t][w][g][a][0] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1];
413 for (
int t=0; t<=1; ++t) {
414 for (
int w=0; w<=1; ++w) {
415 for (
int g=0; g<=1; ++g) {
416 nn[t][w][g][0][0] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0];
423 for (
int t=0; t<=1; ++t) {
424 for (
int w=0; w<=1; ++w) {
425 nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0];
431 for (
int t=0; t<=1; ++t) {
432 nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0];
437 nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0];
439 return nn[0][0][0][0][0];
450 inline VFPEvaluation interpolate(
451 const VFPInjTable::array_type& array,
452 const InterpData& flo_i,
453 const InterpData& thp_i) {
456 VFPEvaluation nn[2][2];
461 for (
int t=0; t<=1; ++t) {
462 for (
int f=0; f<=1; ++f) {
464 const int ti = thp_i.ind_[t];
465 const int fi = flo_i.ind_[f];
468 nn[t][f].value = array[ti][fi];
475 for (
int i=0; i<=1; ++i) {
476 nn[0][i].dthp = (nn[1][i].value - nn[0][i].value) * thp_i.inv_dist_;
477 nn[i][0].dwfr = -1e100;
478 nn[i][0].dgfr = -1e100;
479 nn[i][0].dalq = -1e100;
480 nn[i][0].dflo = (nn[i][1].value - nn[i][0].value) * flo_i.inv_dist_;
482 nn[1][i].dthp = nn[0][i].dthp;
483 nn[i][1].dwfr = nn[i][0].dwfr;
484 nn[i][1].dgfr = nn[i][0].dgfr;
485 nn[i][1].dalq = nn[i][0].dalq;
486 nn[i][1].dflo = nn[i][0].dflo;
494 nn[0][0] = t1*nn[0][0] + t2*nn[0][1];
495 nn[1][0] = t1*nn[1][0] + t2*nn[1][1];
500 nn[0][0] = t1*nn[0][0] + t2*nn[1][0];
509 #pragma GCC pop_options //unroll loops
516 inline VFPEvaluation bhp(
const VFPProdTable* table,
518 const double& liquid,
519 const double& vapour,
523 double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
524 double wfr = detail::getWFR(aqua, liquid, vapour, table->getWFRType());
525 double gfr = detail::getGFR(aqua, liquid, vapour, table->getGFRType());
529 auto flo_i = detail::findInterpData(-flo, table->getFloAxis());
530 auto thp_i = detail::findInterpData( thp, table->getTHPAxis());
531 auto wfr_i = detail::findInterpData( wfr, table->getWFRAxis());
532 auto gfr_i = detail::findInterpData( gfr, table->getGFRAxis());
533 auto alq_i = detail::findInterpData( alq, table->getALQAxis());
535 detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
544 inline VFPEvaluation bhp(
const VFPInjTable* table,
546 const double& liquid,
547 const double& vapour,
550 double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
553 auto flo_i = detail::findInterpData(flo, table->getFloAxis());
554 auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
557 detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i);
572 template <
typename T>
573 const T* getTable(
const std::map<int, T*> tables,
int table_id) {
574 auto entry = tables.find(table_id);
575 if (entry == tables.end()) {
576 OPM_THROW(std::invalid_argument,
"Nonexistent table " << table_id <<
" referenced.");
579 return entry->second;
595 inline void extendBlockPattern(
const ADB& x, std::vector<int>& block_pattern) {
596 std::vector<int> x_block_pattern = x.blockPattern();
598 if (x_block_pattern.empty()) {
602 if (block_pattern.empty()) {
603 block_pattern = x_block_pattern;
607 if (x_block_pattern != block_pattern) {
608 OPM_THROW(std::logic_error,
"Block patterns do not match");
617 inline std::vector<int> commonBlockPattern(
622 std::vector<int> block_pattern;
624 extendBlockPattern(x1, block_pattern);
625 extendBlockPattern(x2, block_pattern);
626 extendBlockPattern(x3, block_pattern);
627 extendBlockPattern(x4, block_pattern);
629 return block_pattern;
632 inline std::vector<int> commonBlockPattern(
638 std::vector<int> block_pattern = commonBlockPattern(x1, x2, x3, x4);
639 extendBlockPattern(x5, block_pattern);
641 return block_pattern;
658 template <
typename TYPE,
typename TABLE>
659 TYPE getType(
const TABLE* table);
663 VFPProdTable::FLO_TYPE getType(
const VFPProdTable* table) {
664 return table->getFloType();
669 VFPProdTable::WFR_TYPE getType(
const VFPProdTable* table) {
670 return table->getWFRType();
675 VFPProdTable::GFR_TYPE getType(
const VFPProdTable* table) {
676 return table->getGFRType();
685 VFPInjTable::FLO_TYPE getType(
const VFPInjTable* table) {
686 return table->getFloType();
695 template <
typename TYPE>
699 const ADB& vapour, TYPE type);
707 VFPProdTable::FLO_TYPE type) {
708 return detail::getFlo(aqua, liquid, vapour, type);
717 VFPProdTable::WFR_TYPE type) {
718 return detail::getWFR(aqua, liquid, vapour, type);
727 VFPProdTable::GFR_TYPE type) {
728 return detail::getGFR(aqua, liquid, vapour, type);
737 VFPInjTable::FLO_TYPE type) {
738 return detail::getFlo(aqua, liquid, vapour, type);
748 template <
typename TYPE,
typename TABLE>
749 ADB combineADBVars(
const std::vector<const TABLE*>& well_tables,
754 const int num_wells =
static_cast<int>(well_tables.size());
755 assert(aqua.size() == num_wells);
756 assert(liquid.size() == num_wells);
757 assert(vapour.size() == num_wells);
760 std::map<TYPE, ADB> map;
763 std::map<TYPE, std::vector<int> > elems;
767 for (
int i=0; i<num_wells; ++i) {
768 const TABLE* table = well_tables[i];
772 TYPE type = getType<TYPE>(table);
776 if (map.find(type) == map.end()) {
777 map.insert(std::pair<TYPE, ADB>(
779 detail::getValue<TYPE>(aqua, liquid, vapour, type)
784 elems[type].push_back(i);
791 for (
const auto& entry : elems) {
792 const auto& key = entry.first;
793 const auto& value = entry.second;
796 assert(map.find(key) != map.end());
797 const ADB& values = map.find(key)->second;
800 const std::vector<int>& current = value;
803 retval = retval +
superset(
subset(values, current), current, values.size());
813 inline double findX(
const double& x0,
818 const double dx = x1 - x0;
819 const double dy = y1 - y0;
831 x = x0 + (y-y0) * (dx/dy);
855 inline double findTHP(
856 const std::vector<double>& bhp_array,
857 const std::vector<double>& thp_array,
859 int nthp = thp_array.size();
864 assert(std::is_sorted(thp_array.begin(), thp_array.end()));
872 if (std::is_sorted(bhp_array.begin(), bhp_array.end())) {
874 if (bhp <= bhp_array[0]) {
876 const double& x0 = thp_array[0];
877 const double& x1 = thp_array[1];
878 const double& y0 = bhp_array[0];
879 const double& y1 = bhp_array[1];
880 thp = detail::findX(x0, x1, y0, y1, bhp);
883 else if (bhp > bhp_array[nthp-1]) {
885 const double& x0 = thp_array[nthp-2];
886 const double& x1 = thp_array[nthp-1];
887 const double& y0 = bhp_array[nthp-2];
888 const double& y1 = bhp_array[nthp-1];
889 thp = detail::findX(x0, x1, y0, y1, bhp);
901 for (; i<nthp-1; ++i) {
902 const double& y0 = bhp_array[i ];
903 const double& y1 = bhp_array[i+1];
905 if (y0 < bhp && bhp <= y1) {
911 assert(found ==
true);
912 static_cast<void>(found);
914 const double& x0 = thp_array[i ];
915 const double& x1 = thp_array[i+1];
916 const double& y0 = bhp_array[i ];
917 const double& y1 = bhp_array[i+1];
918 thp = detail::findX(x0, x1, y0, y1, bhp);
928 for (; i<nthp-1; ++i) {
929 const double& y0 = bhp_array[i ];
930 const double& y1 = bhp_array[i+1];
932 if (y0 < bhp && bhp <= y1) {
938 const double& x0 = thp_array[i ];
939 const double& x1 = thp_array[i+1];
940 const double& y0 = bhp_array[i ];
941 const double& y1 = bhp_array[i+1];
942 thp = detail::findX(x0, x1, y0, y1, bhp);
944 else if (bhp <= bhp_array[0]) {
946 const double& x0 = thp_array[0];
947 const double& x1 = thp_array[1];
948 const double& y0 = bhp_array[0];
949 const double& y1 = bhp_array[1];
950 thp = detail::findX(x0, x1, y0, y1, bhp);
953 else if (bhp > bhp_array[nthp-1]) {
955 const double& x0 = thp_array[nthp-2];
956 const double& x1 = thp_array[nthp-1];
957 const double& y0 = bhp_array[nthp-2];
958 const double& y1 = bhp_array[nthp-1];
959 thp = detail::findX(x0, x1, y0, y1, bhp);
962 OPM_THROW(std::logic_error,
"Programmer error: Unable to find THP in THP array");
An "ADB-like" structure with a single value and a set of derivatives.
Definition: VFPHelpers.hpp:272
Helper struct for linear interpolation.
Definition: VFPHelpers.hpp:192
AutoDiffBlock< Scalar > superset(const AutoDiffBlock< Scalar > &x, const IntVec &indices, const int n)
Returns v where v(indices) == x, v(!indices) == 0 and v.size() == n.
Definition: AutoDiffHelpers.hpp:319
static AutoDiffBlock constant(V &&val)
Create an AutoDiffBlock representing a constant.
Definition: AutoDiffBlock.hpp:111
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:292
Eigen::Array< double, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:99