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>
52 return (std::isnan(value)) ? 0.0 : value;
68 ADB retval = not_nan_selector.select(values, zero);
82 static T
getFlo(
const T& aqua,
const T& liquid,
const T& vapour,
83 const VFPProdTable::FLO_TYPE& type) {
85 case VFPProdTable::FLO_OIL:
88 case VFPProdTable::FLO_LIQ:
91 case VFPProdTable::FLO_GAS:
94 case VFPProdTable::FLO_INVALID:
96 OPM_THROW(std::logic_error,
"Invalid FLO_TYPE: '" << type <<
"'");
107 template <
typename T>
108 static T
getFlo(
const T& aqua,
const T& liquid,
const T& vapour,
109 const VFPInjTable::FLO_TYPE& type) {
111 case VFPInjTable::FLO_OIL:
114 case VFPInjTable::FLO_WAT:
117 case VFPInjTable::FLO_GAS:
120 case VFPInjTable::FLO_INVALID:
122 OPM_THROW(std::logic_error,
"Invalid FLO_TYPE: '" << type <<
"'");
136 template <
typename T>
137 static T
getWFR(
const T& aqua,
const T& liquid,
const T& vapour,
138 const VFPProdTable::WFR_TYPE& type) {
140 case VFPProdTable::WFR_WOR: {
142 T wor = aqua / liquid;
145 case VFPProdTable::WFR_WCT:
147 return zeroIfNan(aqua / (aqua + liquid));
148 case VFPProdTable::WFR_WGR:
151 case VFPProdTable::WFR_INVALID:
153 OPM_THROW(std::logic_error,
"Invalid WFR_TYPE: '" << type <<
"'");
166 template <
typename T>
167 static T
getGFR(
const T& aqua,
const T& liquid,
const T& vapour,
168 const VFPProdTable::GFR_TYPE& type) {
170 case VFPProdTable::GFR_GOR:
173 case VFPProdTable::GFR_GLR:
175 return zeroIfNan(vapour / (liquid + aqua));
176 case VFPProdTable::GFR_OGR:
179 case VFPProdTable::GFR_INVALID:
181 OPM_THROW(std::logic_error,
"Invalid GFR_TYPE: '" << type <<
"'");
214 const int nvalues = values.size();
217 if (values.size() == 1) {
226 if (value < values.front()) {
231 else if (value >= values.back()) {
232 retval.
ind_[0] = nvalues-2;
233 retval.
ind_[1] = nvalues-1;
237 for (
int i=1; i<nvalues; ++i) {
238 if (values[i] >= value) {
239 retval.
ind_[0] = i-1;
246 const double start = values[retval.
ind_[0]];
247 const double end = values[retval.
ind_[1]];
329 #pragma GCC push_options
330 #pragma GCC optimize ("unroll-loops")
335 const VFPProdTable::array_type& array,
351 for (
int t=0; t<=1; ++t) {
352 for (
int w=0; w<=1; ++w) {
353 for (
int g=0; g<=1; ++g) {
354 for (
int a=0; a<=1; ++a) {
355 for (
int f=0; f<=1; ++f) {
357 const int ti = thp_i.
ind_[t];
358 const int wi = wfr_i.
ind_[w];
359 const int gi = gfr_i.
ind_[g];
360 const int ai = alq_i.
ind_[a];
361 const int fi = flo_i.
ind_[f];
364 nn[t][w][g][a][f].
value = array[ti][wi][gi][ai][fi];
374 for (
int i=0; i<=1; ++i) {
375 for (
int j=0; j<=1; ++j) {
376 for (
int k=0; k<=1; ++k) {
377 for (
int l=0; l<=1; ++l) {
384 nn[1][i][j][k][l].
dthp = nn[0][i][j][k][l].
dthp;
385 nn[i][1][j][k][l].
dwfr = nn[i][0][j][k][l].
dwfr;
386 nn[i][j][1][k][l].
dgfr = nn[i][j][0][k][l].
dgfr;
387 nn[i][j][k][1][l].
dalq = nn[i][j][k][0][l].
dalq;
388 nn[i][j][k][l][1].
dflo = nn[i][j][k][l][0].
dflo;
402 for (
int t=0; t<=1; ++t) {
403 for (
int w=0; w<=1; ++w) {
404 for (
int g=0; g<=1; ++g) {
405 for (
int a=0; a<=1; ++a) {
406 nn[t][w][g][a][0] = t1*nn[t][w][g][a][0] + t2*nn[t][w][g][a][1];
414 for (
int t=0; t<=1; ++t) {
415 for (
int w=0; w<=1; ++w) {
416 for (
int g=0; g<=1; ++g) {
417 nn[t][w][g][0][0] = t1*nn[t][w][g][0][0] + t2*nn[t][w][g][1][0];
424 for (
int t=0; t<=1; ++t) {
425 for (
int w=0; w<=1; ++w) {
426 nn[t][w][0][0][0] = t1*nn[t][w][0][0][0] + t2*nn[t][w][1][0][0];
432 for (
int t=0; t<=1; ++t) {
433 nn[t][0][0][0][0] = t1*nn[t][0][0][0][0] + t2*nn[t][1][0][0][0];
438 nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0];
440 return nn[0][0][0][0][0];
452 const VFPInjTable::array_type& array,
462 for (
int t=0; t<=1; ++t) {
463 for (
int f=0; f<=1; ++f) {
465 const int ti = thp_i.
ind_[t];
466 const int fi = flo_i.
ind_[f];
469 nn[t][f].
value = array[ti][fi];
476 for (
int i=0; i<=1; ++i) {
478 nn[i][0].
dwfr = -1e100;
479 nn[i][0].
dgfr = -1e100;
480 nn[i][0].
dalq = -1e100;
495 nn[0][0] = t1*nn[0][0] + t2*nn[0][1];
496 nn[1][0] = t1*nn[1][0] + t2*nn[1][1];
501 nn[0][0] = t1*nn[0][0] + t2*nn[1][0];
510 #pragma GCC pop_options //unroll loops
519 const double& liquid,
520 const double& vapour,
524 double flo =
detail::getFlo(aqua, liquid, vapour, table->getFloType());
525 double wfr =
detail::getWFR(aqua, liquid, vapour, table->getWFRType());
526 double gfr =
detail::getGFR(aqua, liquid, vapour, table->getGFRType());
547 const double& liquid,
548 const double& vapour,
551 double flo =
detail::getFlo(aqua, liquid, vapour, table->getFloType());
573 template <
typename T>
574 const T*
getTable(
const std::map<int, T*> tables,
int table_id) {
575 auto entry = tables.find(table_id);
576 if (entry == tables.end()) {
577 OPM_THROW(std::invalid_argument,
"Nonexistent table " << table_id <<
" referenced.");
580 return entry->second;
599 if (x_block_pattern.empty()) {
603 if (block_pattern.empty()) {
604 block_pattern = x_block_pattern;
608 if (x_block_pattern != block_pattern) {
609 OPM_THROW(std::logic_error,
"Block patterns do not match");
623 std::vector<int> block_pattern;
630 return block_pattern;
642 return block_pattern;
659 template <
typename TYPE,
typename TABLE>
660 TYPE
getType(
const TABLE* table);
664 VFPProdTable::FLO_TYPE
getType(
const VFPProdTable* table) {
665 return table->getFloType();
670 VFPProdTable::WFR_TYPE
getType(
const VFPProdTable* table) {
671 return table->getWFRType();
676 VFPProdTable::GFR_TYPE
getType(
const VFPProdTable* table) {
677 return table->getGFRType();
686 VFPInjTable::FLO_TYPE
getType(
const VFPInjTable* table) {
687 return table->getFloType();
696 template <
typename TYPE>
700 const ADB& vapour, TYPE type);
708 VFPProdTable::FLO_TYPE type) {
718 VFPProdTable::WFR_TYPE type) {
728 VFPProdTable::GFR_TYPE type) {
738 VFPInjTable::FLO_TYPE type) {
749 template <
typename TYPE,
typename TABLE>
755 const int num_wells =
static_cast<int>(well_tables.size());
756 assert(aqua.
size() == num_wells);
757 assert(liquid.
size() == num_wells);
758 assert(vapour.
size() == num_wells);
761 std::map<TYPE, ADB> map;
764 std::map<TYPE, std::vector<int> > elems;
768 for (
int i=0; i<num_wells; ++i) {
769 const TABLE* table = well_tables[i];
773 TYPE type = getType<TYPE>(table);
777 if (map.find(type) == map.end()) {
778 map.insert(std::pair<TYPE, ADB>(
780 detail::getValue<TYPE>(aqua, liquid, vapour, type)
785 elems[type].push_back(i);
792 for (
const auto& entry : elems) {
793 const auto& key = entry.first;
794 const auto& value = entry.second;
797 assert(map.find(key) != map.end());
798 const ADB& values = map.find(key)->second;
801 const std::vector<int>& current = value;
804 retval = retval +
superset(
subset(values, current), current, values.size());
814 inline double findX(
const double& x0,
819 const double dx = x1 - x0;
820 const double dy = y1 - y0;
832 x = x0 + (y-y0) * (dx/dy);
857 const std::vector<double>& bhp_array,
858 const std::vector<double>& thp_array,
860 int nthp = thp_array.size();
865 assert(std::is_sorted(thp_array.begin(), thp_array.end()));
873 if (std::is_sorted(bhp_array.begin(), bhp_array.end())) {
875 if (bhp <= bhp_array[0]) {
877 const double& x0 = thp_array[0];
878 const double& x1 = thp_array[1];
879 const double& y0 = bhp_array[0];
880 const double& y1 = bhp_array[1];
884 else if (bhp > bhp_array[nthp-1]) {
886 const double& x0 = thp_array[nthp-2];
887 const double& x1 = thp_array[nthp-1];
888 const double& y0 = bhp_array[nthp-2];
889 const double& y1 = bhp_array[nthp-1];
902 for (; i<nthp-1; ++i) {
903 const double& y0 = bhp_array[i ];
904 const double& y1 = bhp_array[i+1];
906 if (y0 < bhp && bhp <= y1) {
912 assert(found ==
true);
913 static_cast<void>(found);
915 const double& x0 = thp_array[i ];
916 const double& x1 = thp_array[i+1];
917 const double& y0 = bhp_array[i ];
918 const double& y1 = bhp_array[i+1];
929 for (; i<nthp-1; ++i) {
930 const double& y0 = bhp_array[i ];
931 const double& y1 = bhp_array[i+1];
933 if (y0 < bhp && bhp <= y1) {
939 const double& x0 = thp_array[i ];
940 const double& x1 = thp_array[i+1];
941 const double& y0 = bhp_array[i ];
942 const double& y1 = bhp_array[i+1];
945 else if (bhp <= bhp_array[0]) {
947 const double& x0 = thp_array[0];
948 const double& x1 = thp_array[1];
949 const double& y0 = bhp_array[0];
950 const double& y1 = bhp_array[1];
954 else if (bhp > bhp_array[nthp-1]) {
956 const double& x0 = thp_array[nthp-2];
957 const double& x1 = thp_array[nthp-1];
958 const double& y0 = bhp_array[nthp-2];
959 const double& y1 = bhp_array[nthp-1];
963 OPM_THROW(std::logic_error,
"Programmer error: Unable to find THP in THP array");
VFPEvaluation bhp(const VFPProdTable *table, const double &aqua, const double &liquid, const double &vapour, const double &thp, const double &alq)
Definition: VFPHelpers.hpp:517
static T getWFR(const T &aqua, const T &liquid, const T &vapour, const VFPProdTable::WFR_TYPE &type)
Definition: VFPHelpers.hpp:137
Eigen::Array< Scalar, Eigen::Dynamic, 1 > subset(const Eigen::Array< Scalar, Eigen::Dynamic, 1 > &x, const IntVec &indices)
Returns x(indices).
Definition: AutoDiffHelpers.hpp:287
double inv_dist_
Definition: VFPHelpers.hpp:196
Eigen::Array< Scalar, Eigen::Dynamic, 1 > V
Underlying type for values.
Definition: AutoDiffBlock.hpp:98
Definition: AdditionalObjectDeleter.hpp:22
std::vector< int > blockPattern() const
Sizes (number of columns) of Jacobian blocks.
Definition: AutoDiffBlock.hpp:425
VFPEvaluation interpolate(const VFPProdTable::array_type &array, const InterpData &flo_i, const InterpData &thp_i, const InterpData &wfr_i, const InterpData &gfr_i, const InterpData &alq_i)
Definition: VFPHelpers.hpp:334
int ind_[2]
Definition: VFPHelpers.hpp:195
std::vector< int > commonBlockPattern(const ADB &x1, const ADB &x2, const ADB &x3, const ADB &x4)
Definition: VFPHelpers.hpp:618
double zeroIfNan(const double &value)
Definition: VFPHelpers.hpp:51
static T getGFR(const T &aqua, const T &liquid, const T &vapour, const VFPProdTable::GFR_TYPE &type)
Definition: VFPHelpers.hpp:167
Selection. Choose first of two elements if selection basis element is nonnegative.
Definition: AutoDiffHelpers.hpp:359
InterpData findInterpData(const double &value, const std::vector< double > &values)
Definition: VFPHelpers.hpp:211
double dthp
Definition: VFPHelpers.hpp:276
int size() const
Number of elements.
Definition: AutoDiffBlock.hpp:413
TYPE getType(const TABLE *table)
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:314
static T getFlo(const T &aqua, const T &liquid, const T &vapour, const VFPProdTable::FLO_TYPE &type)
Definition: VFPHelpers.hpp:82
double findTHP(const std::vector< double > &bhp_array, const std::vector< double > &thp_array, double bhp)
Definition: VFPHelpers.hpp:856
static AutoDiffBlock constant(V &&val)
Definition: AutoDiffBlock.hpp:110
InterpData()
Definition: VFPHelpers.hpp:194
double dalq
Definition: VFPHelpers.hpp:279
ADB combineADBVars(const std::vector< const TABLE * > &well_tables, const ADB &aqua, const ADB &liquid, const ADB &vapour)
Definition: VFPHelpers.hpp:750
Definition: VFPHelpers.hpp:273
void extendBlockPattern(const ADB &x, std::vector< int > &block_pattern)
Definition: VFPHelpers.hpp:596
double dgfr
Definition: VFPHelpers.hpp:278
double findX(const double &x0, const double &x1, const double &y0, const double &y1, const double &y)
Definition: VFPHelpers.hpp:814
double value
Definition: VFPHelpers.hpp:274
VFPEvaluation operator*(double lhs, const VFPEvaluation &rhs)
Definition: VFPHelpers.hpp:307
ADB getValue(const ADB &aqua, const ADB &liquid, const ADB &vapour, TYPE type)
Definition: VFPHelpers.hpp:193
VFPEvaluation operator-(VFPEvaluation lhs, const VFPEvaluation &rhs)
Definition: VFPHelpers.hpp:295
const V & value() const
Function value.
Definition: AutoDiffBlock.hpp:436
const T * getTable(const std::map< int, T * > tables, int table_id)
Definition: VFPHelpers.hpp:574
AutoDiffBlock< double > ADB
Definition: VFPHelpers.hpp:37
double factor_
Definition: VFPHelpers.hpp:197
VFPEvaluation operator+(VFPEvaluation lhs, const VFPEvaluation &rhs)
Definition: VFPHelpers.hpp:283
VFPEvaluation()
Definition: VFPHelpers.hpp:274
double dwfr
Definition: VFPHelpers.hpp:277
double dflo
Definition: VFPHelpers.hpp:280