VFPHelpers.hpp
1 /*
2  Copyright 2015 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 
21 #ifndef OPM_AUTODIFF_VFPHELPERS_HPP_
22 #define OPM_AUTODIFF_VFPHELPERS_HPP_
23 
24 
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>
30 
34 namespace Opm {
35 namespace detail {
36 
37 
38 typedef AutoDiffBlock<double> ADB;
39 
40 
44 inline double zeroIfNan(const double& value) {
45  return (std::isnan(value)) ? 0.0 : value;
46 }
47 
51 template <class EvalWell>
52 inline double zeroIfNan(const EvalWell& value) {
53  return (std::isnan(value.value())) ? 0.0 : value.value();
54 }
55 
56 
57 
61 inline ADB zeroIfNan(const ADB& values) {
62  Selector<ADB::V::Scalar> not_nan_selector(values.value(), Selector<ADB::V::Scalar>::NotNaN);
63 
64  const ADB::V z = ADB::V::Zero(values.size());
65  const ADB zero = ADB::constant(z, values.blockPattern());
66 
67  ADB retval = not_nan_selector.select(values, zero);
68  return retval;
69 }
70 
71 
72 
73 
74 
80 template <typename T>
81 static T getFlo(const T& aqua, const T& liquid, const T& vapour,
82  const VFPProdTable::FLO_TYPE& type) {
83  switch (type) {
84  case VFPProdTable::FLO_OIL:
85  //Oil = liquid phase
86  return liquid;
87  case VFPProdTable::FLO_LIQ:
88  //Liquid = aqua + liquid phases
89  return aqua + liquid;
90  case VFPProdTable::FLO_GAS:
91  //Gas = vapor phase
92  return vapour;
93  case VFPProdTable::FLO_INVALID: //Intentional fall-through
94  default:
95  OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
96  }
97 }
98 
99 
100 
106 template <typename T>
107 static T getFlo(const T& aqua, const T& liquid, const T& vapour,
108  const VFPInjTable::FLO_TYPE& type) {
109  switch (type) {
110  case VFPInjTable::FLO_OIL:
111  //Oil = liquid phase
112  return liquid;
113  case VFPInjTable::FLO_WAT:
114  //Liquid = aqua phase
115  return aqua;
116  case VFPInjTable::FLO_GAS:
117  //Gas = vapor phase
118  return vapour;
119  case VFPInjTable::FLO_INVALID: //Intentional fall-through
120  default:
121  OPM_THROW(std::logic_error, "Invalid FLO_TYPE: '" << type << "'");
122  }
123 }
124 
125 
126 
127 
128 
129 
130 
135 template <typename T>
136 static T getWFR(const T& aqua, const T& liquid, const T& vapour,
137  const VFPProdTable::WFR_TYPE& type) {
138  switch(type) {
139  case VFPProdTable::WFR_WOR: {
140  //Water-oil ratio = water / oil
141  T wor = aqua / liquid;
142  return zeroIfNan(wor);
143  }
144  case VFPProdTable::WFR_WCT:
145  //Water cut = water / (water + oil)
146  return zeroIfNan(aqua / (aqua + liquid));
147  case VFPProdTable::WFR_WGR:
148  //Water-gas ratio = water / gas
149  return zeroIfNan(aqua / vapour);
150  case VFPProdTable::WFR_INVALID: //Intentional fall-through
151  default:
152  OPM_THROW(std::logic_error, "Invalid WFR_TYPE: '" << type << "'");
153  }
154 }
155 
156 
157 
158 
159 
160 
165 template <typename T>
166 static T getGFR(const T& aqua, const T& liquid, const T& vapour,
167  const VFPProdTable::GFR_TYPE& type) {
168  switch(type) {
169  case VFPProdTable::GFR_GOR:
170  // Gas-oil ratio = gas / oil
171  return zeroIfNan(vapour / liquid);
172  case VFPProdTable::GFR_GLR:
173  // Gas-liquid ratio = gas / (oil + water)
174  return zeroIfNan(vapour / (liquid + aqua));
175  case VFPProdTable::GFR_OGR:
176  // Oil-gas ratio = oil / gas
177  return zeroIfNan(liquid / vapour);
178  case VFPProdTable::GFR_INVALID: //Intentional fall-through
179  default:
180  OPM_THROW(std::logic_error, "Invalid GFR_TYPE: '" << type << "'");
181  }
182 }
183 
184 
185 
186 
187 
188 
192 struct InterpData {
193  InterpData() : ind_{0, 0}, inv_dist_(0.0), factor_(0.0) {}
194  int ind_[2]; //[First element greater than or equal to value, Last element smaller than or equal to value]
195  double inv_dist_; // 1 / distance between the two end points of the segment. Used to calculate derivatives and uses 1.0 / 0.0 = 0.0 as a convention
196  double factor_; // Interpolation factor
197 };
198 
199 
200 
201 
202 
203 
210 inline InterpData findInterpData(const double& value, const std::vector<double>& values) {
211  InterpData retval;
212 
213  const int nvalues = values.size();
214 
215  //If we only have one value in our vector, return that
216  if (values.size() == 1) {
217  retval.ind_[0] = 0;
218  retval.ind_[1] = 0;
219  retval.inv_dist_ = 0.0;
220  retval.factor_ = 0.0;
221  }
222  // Else search in the vector
223  else {
224  //If value is less than all values, use first interval
225  if (value < values.front()) {
226  retval.ind_[0] = 0;
227  retval.ind_[1] = 1;
228  }
229  //If value is greater than all values, use last interval
230  else if (value >= values.back()) {
231  retval.ind_[0] = nvalues-2;
232  retval.ind_[1] = nvalues-1;
233  }
234  else {
235  //Search internal intervals
236  for (int i=1; i<nvalues; ++i) {
237  if (values[i] >= value) {
238  retval.ind_[0] = i-1;
239  retval.ind_[1] = i;
240  break;
241  }
242  }
243  }
244 
245  const double start = values[retval.ind_[0]];
246  const double end = values[retval.ind_[1]];
247 
248  //Find interpolation ratio
249  if (end > start) {
250  //FIXME: Possible source for floating point error here if value and floor are large,
251  //but very close to each other
252  retval.inv_dist_ = 1.0 / (end-start);
253  retval.factor_ = (value-start) * retval.inv_dist_;
254  }
255  else {
256  retval.inv_dist_ = 0.0;
257  retval.factor_ = 0.0;
258  }
259  }
260 
261  return retval;
262 }
263 
264 
265 
266 
267 
268 
273  VFPEvaluation() : value(0.0), dthp(0.0), dwfr(0.0), dgfr(0.0), dalq(0.0), dflo(0.0) {};
274  double value;
275  double dthp;
276  double dwfr;
277  double dgfr;
278  double dalq;
279  double dflo;
280 };
281 
282 inline VFPEvaluation operator+(
283  VFPEvaluation lhs,
284  const VFPEvaluation& rhs) {
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;
291  return lhs;
292 }
293 
294 inline VFPEvaluation operator-(
295  VFPEvaluation lhs,
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;
303  return lhs;
304 }
305 
306 inline VFPEvaluation operator*(
307  double lhs,
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;
316  return retval;
317 }
318 
319 
320 
321 
322 
323 
327 #ifdef __GNUC__
328 #pragma GCC push_options
329 #pragma GCC optimize ("unroll-loops")
330 #endif
331 
332 
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) {
340 
341  //Values and derivatives in a 5D hypercube
342  VFPEvaluation nn[2][2][2][2][2];
343 
344 
345  //Pick out nearest neighbors (nn) to our evaluation point
346  //This is not really required, but performance-wise it may pay off, since the 32-elements
347  //we copy to (nn) will fit better in cache than the full original table for the
348  //interpolation below.
349  //The following ladder of for loops will presumably be unrolled by a reasonable compiler.
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) {
355  //Shorthands for indexing
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];
361 
362  //Copy element
363  nn[t][w][g][a][f].value = array[ti][wi][gi][ai][fi];
364  }
365  }
366  }
367  }
368  }
369 
370  //Calculate derivatives
371  //Note that the derivative of the two end points of a line aligned with the
372  //"axis of the derivative" are equal
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_;
382 
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;
388  }
389  }
390  }
391  }
392 
393  double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
394 
395  // Remove dimensions one by one
396  // Example: going from 3D to 2D to 1D, we start by interpolating along
397  // the z axis first, leaving a 2D problem. Then interpolating along the y
398  // axis, leaving a 1D, problem, etc.
399  t2 = flo_i.factor_;
400  t1 = (1.0-t2);
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];
406  }
407  }
408  }
409  }
410 
411  t2 = alq_i.factor_;
412  t1 = (1.0-t2);
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];
417  }
418  }
419  }
420 
421  t2 = gfr_i.factor_;
422  t1 = (1.0-t2);
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];
426  }
427  }
428 
429  t2 = wfr_i.factor_;
430  t1 = (1.0-t2);
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];
433  }
434 
435  t2 = thp_i.factor_;
436  t1 = (1.0-t2);
437  nn[0][0][0][0][0] = t1*nn[0][0][0][0][0] + t2*nn[1][0][0][0][0];
438 
439  return nn[0][0][0][0][0];
440 }
441 
442 
443 
444 
445 
450 inline VFPEvaluation interpolate(
451  const VFPInjTable::array_type& array,
452  const InterpData& flo_i,
453  const InterpData& thp_i) {
454 
455  //Values and derivatives in a 2D plane
456  VFPEvaluation nn[2][2];
457 
458 
459  //Pick out nearest neighbors (nn) to our evaluation point
460  //The following ladder of for loops will presumably be unrolled by a reasonable compiler.
461  for (int t=0; t<=1; ++t) {
462  for (int f=0; f<=1; ++f) {
463  //Shorthands for indexing
464  const int ti = thp_i.ind_[t];
465  const int fi = flo_i.ind_[f];
466 
467  //Copy element
468  nn[t][f].value = array[ti][fi];
469  }
470  }
471 
472  //Calculate derivatives
473  //Note that the derivative of the two end points of a line aligned with the
474  //"axis of the derivative" are equal
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_;
481 
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;
487  }
488 
489  double t1, t2; //interpolation variables, so that t1 = (1-t) and t2 = t.
490 
491  // Go from 2D to 1D
492  t2 = flo_i.factor_;
493  t1 = (1.0-t2);
494  nn[0][0] = t1*nn[0][0] + t2*nn[0][1];
495  nn[1][0] = t1*nn[1][0] + t2*nn[1][1];
496 
497  // Go from line to point on line
498  t2 = thp_i.factor_;
499  t1 = (1.0-t2);
500  nn[0][0] = t1*nn[0][0] + t2*nn[1][0];
501 
502  return nn[0][0];
503 }
504 
505 
506 
507 
508 #ifdef __GNUC__
509 #pragma GCC pop_options //unroll loops
510 #endif
511 
512 
513 
514 
515 
516 inline VFPEvaluation bhp(const VFPProdTable* table,
517  const double& aqua,
518  const double& liquid,
519  const double& vapour,
520  const double& thp,
521  const double& alq) {
522  //Find interpolation variables
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());
526 
527  //First, find the values to interpolate between
528  //Recall that flo is negative in Opm, so switch sign.
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());
534 
535  detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i, wfr_i, gfr_i, alq_i);
536 
537  return retval;
538 }
539 
540 
541 
542 
543 
544 inline VFPEvaluation bhp(const VFPInjTable* table,
545  const double& aqua,
546  const double& liquid,
547  const double& vapour,
548  const double& thp) {
549  //Find interpolation variables
550  double flo = detail::getFlo(aqua, liquid, vapour, table->getFloType());
551 
552  //First, find the values to interpolate between
553  auto flo_i = detail::findInterpData(flo, table->getFloAxis());
554  auto thp_i = detail::findInterpData(thp, table->getTHPAxis());
555 
556  //Then perform the interpolation itself
557  detail::VFPEvaluation retval = detail::interpolate(table->getTable(), flo_i, thp_i);
558 
559  return retval;
560 }
561 
562 
563 
564 
565 
566 
567 
568 
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.");
577  }
578  else {
579  return entry->second;
580  }
581 }
582 
583 
584 
585 
586 
587 
588 
589 
590 
591 
595 inline void extendBlockPattern(const ADB& x, std::vector<int>& block_pattern) {
596  std::vector<int> x_block_pattern = x.blockPattern();
597 
598  if (x_block_pattern.empty()) {
599  return;
600  }
601  else {
602  if (block_pattern.empty()) {
603  block_pattern = x_block_pattern;
604  return;
605  }
606  else {
607  if (x_block_pattern != block_pattern) {
608  OPM_THROW(std::logic_error, "Block patterns do not match");
609  }
610  }
611  }
612 }
613 
617 inline std::vector<int> commonBlockPattern(
618  const ADB& x1,
619  const ADB& x2,
620  const ADB& x3,
621  const ADB& x4) {
622  std::vector<int> block_pattern;
623 
624  extendBlockPattern(x1, block_pattern);
625  extendBlockPattern(x2, block_pattern);
626  extendBlockPattern(x3, block_pattern);
627  extendBlockPattern(x4, block_pattern);
628 
629  return block_pattern;
630 }
631 
632 inline std::vector<int> commonBlockPattern(
633  const ADB& x1,
634  const ADB& x2,
635  const ADB& x3,
636  const ADB& x4,
637  const ADB& x5) {
638  std::vector<int> block_pattern = commonBlockPattern(x1, x2, x3, x4);
639  extendBlockPattern(x5, block_pattern);
640 
641  return block_pattern;
642 }
643 
644 
645 
646 
647 
648 
649 
650 
651 
652 
653 
654 
658 template <typename TYPE, typename TABLE>
659 TYPE getType(const TABLE* table);
660 
661 template <>
662 inline
663 VFPProdTable::FLO_TYPE getType(const VFPProdTable* table) {
664  return table->getFloType();
665 }
666 
667 template <>
668 inline
669 VFPProdTable::WFR_TYPE getType(const VFPProdTable* table) {
670  return table->getWFRType();
671 }
672 
673 template <>
674 inline
675 VFPProdTable::GFR_TYPE getType(const VFPProdTable* table) {
676  return table->getGFRType();
677 }
678 
679 
683 template <>
684 inline
685 VFPInjTable::FLO_TYPE getType(const VFPInjTable* table) {
686  return table->getFloType();
687 }
688 
689 
690 
691 
695 template <typename TYPE>
696 ADB getValue(
697  const ADB& aqua,
698  const ADB& liquid,
699  const ADB& vapour, TYPE type);
700 
701 template <>
702 inline
703 ADB getValue(
704  const ADB& aqua,
705  const ADB& liquid,
706  const ADB& vapour,
707  VFPProdTable::FLO_TYPE type) {
708  return detail::getFlo(aqua, liquid, vapour, type);
709 }
710 
711 template <>
712 inline
713 ADB getValue(
714  const ADB& aqua,
715  const ADB& liquid,
716  const ADB& vapour,
717  VFPProdTable::WFR_TYPE type) {
718  return detail::getWFR(aqua, liquid, vapour, type);
719 }
720 
721 template <>
722 inline
723 ADB getValue(
724  const ADB& aqua,
725  const ADB& liquid,
726  const ADB& vapour,
727  VFPProdTable::GFR_TYPE type) {
728  return detail::getGFR(aqua, liquid, vapour, type);
729 }
730 
731 template <>
732 inline
733 ADB getValue(
734  const ADB& aqua,
735  const ADB& liquid,
736  const ADB& vapour,
737  VFPInjTable::FLO_TYPE type) {
738  return detail::getFlo(aqua, liquid, vapour, type);
739 }
740 
748 template <typename TYPE, typename TABLE>
749 ADB combineADBVars(const std::vector<const TABLE*>& well_tables,
750  const ADB& aqua,
751  const ADB& liquid,
752  const ADB& vapour) {
753 
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);
758 
759  //Caching variable for flo/wfr/gfr
760  std::map<TYPE, ADB> map;
761 
762  //Indexing variable used when combining the different ADB types
763  std::map<TYPE, std::vector<int> > elems;
764 
765  //Compute all of the different ADB types,
766  //and record which wells use which types
767  for (int i=0; i<num_wells; ++i) {
768  const TABLE* table = well_tables[i];
769 
770  //Only do something if this well is under THP control
771  if (table != NULL) {
772  TYPE type = getType<TYPE>(table);
773 
774  //"Caching" of flo_type etc: Only calculate used types
775  //Create type if it does not exist
776  if (map.find(type) == map.end()) {
777  map.insert(std::pair<TYPE, ADB>(
778  type,
779  detail::getValue<TYPE>(aqua, liquid, vapour, type)
780  ));
781  }
782 
783  //Add the index for assembly later in gather_vars
784  elems[type].push_back(i);
785  }
786  }
787 
788  //Loop over all types of ADB variables, and combine them
789  //so that each well gets the proper variable
790  ADB retval = ADB::constant(ADB::V::Zero(num_wells));
791  for (const auto& entry : elems) {
792  const auto& key = entry.first;
793  const auto& value = entry.second;
794 
795  //Get the ADB for this type of variable
796  assert(map.find(key) != map.end());
797  const ADB& values = map.find(key)->second;
798 
799  //Get indices to all elements that should use this ADB
800  const std::vector<int>& current = value;
801 
802  //Add these elements to retval
803  retval = retval + superset(subset(values, current), current, values.size());
804  }
805 
806  return retval;
807 }
808 
813 inline double findX(const double& x0,
814  const double& x1,
815  const double& y0,
816  const double& y1,
817  const double& y) {
818  const double dx = x1 - x0;
819  const double dy = y1 - y0;
820 
828  double x = 0.0;
829 
830  if (dy != 0.0) {
831  x = x0 + (y-y0) * (dx/dy);
832  }
833  else {
834  x = x1;
835  }
836 
837  return x;
838 }
839 
840 
841 
842 
843 
844 
845 
846 
847 
848 
855 inline double findTHP(
856  const std::vector<double>& bhp_array,
857  const std::vector<double>& thp_array,
858  double bhp) {
859  int nthp = thp_array.size();
860 
861  double thp = -1e100;
862 
863  //Check that our thp axis is sorted
864  assert(std::is_sorted(thp_array.begin(), thp_array.end()));
865 
872  if (std::is_sorted(bhp_array.begin(), bhp_array.end())) {
873  //Target bhp less than all values in array, extrapolate
874  if (bhp <= bhp_array[0]) {
875  //TODO: LOG extrapolation
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);
881  }
882  //Target bhp greater than all values in array, extrapolate
883  else if (bhp > bhp_array[nthp-1]) {
884  //TODO: LOG extrapolation
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);
890  }
891  //Target bhp within table ranges, interpolate
892  else {
893  //Loop over the values and find min(bhp_array(thp)) == bhp
894  //so that we maximize the rate.
895 
896  //Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
897  //Assuming a small number of values in bhp_array, this should be quite
898  //efficient. Other strategies might be bisection, etc.
899  int i=0;
900  bool found = false;
901  for (; i<nthp-1; ++i) {
902  const double& y0 = bhp_array[i ];
903  const double& y1 = bhp_array[i+1];
904 
905  if (y0 < bhp && bhp <= y1) {
906  found = true;
907  break;
908  }
909  }
910  //Canary in a coal mine: shouldn't really be required
911  assert(found == true);
912  static_cast<void>(found); //Silence compiler warning
913 
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);
919  }
920  }
921  //bhp_array not sorted, raw search.
922  else {
923  //Find i so that bhp_array[i-1] <= bhp <= bhp_array[i];
924  //Since the BHP values might not be sorted, first search within
925  //our interpolation values, and then try to extrapolate.
926  int i=0;
927  bool found = false;
928  for (; i<nthp-1; ++i) {
929  const double& y0 = bhp_array[i ];
930  const double& y1 = bhp_array[i+1];
931 
932  if (y0 < bhp && bhp <= y1) {
933  found = true;
934  break;
935  }
936  }
937  if (found) {
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);
943  }
944  else if (bhp <= bhp_array[0]) {
945  //TODO: LOG extrapolation
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);
951  }
952  //Target bhp greater than all values in array, extrapolate
953  else if (bhp > bhp_array[nthp-1]) {
954  //TODO: LOG extrapolation
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);
960  }
961  else {
962  OPM_THROW(std::logic_error, "Programmer error: Unable to find THP in THP array");
963  }
964  }
965 
966  return thp;
967 }
968 
969 
970 
971 
972 
973 
974 
975 } // namespace detail
976 
977 
978 } // namespace
979 
980 
981 
982 
983 #endif /* OPM_AUTODIFF_VFPHELPERS_HPP_ */
An "ADB-like" structure with a single value and a set of derivatives.
Definition: VFPHelpers.hpp:272
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition: AdditionalObjectDeleter.hpp:22
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