EquilibrationHelpers.hpp
Go to the documentation of this file.
1 /*
2  Copyright 2014 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 #ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
21 #define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
22 
28 
29 #include <memory>
30 
31 
32 /*
33 ---- synopsis of EquilibrationHelpers.hpp ----
34 
35 namespace Opm
36 {
37  namespace Equil {
38 
39  template <class Props>
40  class DensityCalculator;
41 
42  template <>
43  class DensityCalculator< BlackoilPropertiesInterface >;
44 
45  namespace Miscibility {
46  class RsFunction;
47  class NoMixing;
48  class RsVD;
49  class RsSatAtContact;
50  }
51 
52  struct EquilRecord;
53 
54  template <class DensCalc>
55  class EquilReg;
56 
57 
58  struct PcEq;
59 
60  inline double satFromPc(const BlackoilPropertiesInterface& props,
61  const int phase,
62  const int cell,
63  const double target_pc,
64  const bool increasing = false);
65  struct PcEqSum
66  inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
67  const int phase1,
68  const int phase2,
69  const int cell,
70  const double target_pc);
71  } // namespace Equil
72 } // namespace Opm
73 
74 ---- end of synopsis of EquilibrationHelpers.hpp ----
75 */
76 
77 
78 namespace Opm
79 {
87  namespace Equil {
88 
89 
90  template <class Props>
92 
100  template <>
102  public:
113  const int c)
114  : props_(props)
115  , c_(1, c)
116  {
117  }
118 
131  std::vector<double>
132  operator()(const double p,
133  const double T,
134  const std::vector<double>& z) const
135  {
136  const int np = props_.numPhases();
137  std::vector<double> A(np * np, 0);
138 
139  assert (z.size() == std::vector<double>::size_type(np));
140 
141  double* dAdp = 0;
142  props_.matrix(1, &p, &T, &z[0], &c_[0], &A[0], dAdp);
143 
144  std::vector<double> rho(np, 0.0);
145  props_.density(1, &A[0], &c_[0], &rho[0]);
146 
147  return rho;
148  }
149 
150  private:
151  const BlackoilPropertiesInterface& props_;
152  const std::vector<int> c_;
153  };
154 
155 
160  namespace Miscibility {
161 
166  {
167  public:
183  virtual double operator()(const double depth,
184  const double press,
185  const double temp,
186  const double sat = 0.0) const = 0;
187  };
188 
189 
193  class NoMixing : public RsFunction {
194  public:
211  double
212  operator()(const double /* depth */,
213  const double /* press */,
214  const double /* temp */,
215  const double /* sat */ = 0.0) const
216  {
217  return 0.0;
218  }
219  };
220 
221 
227  class RsVD : public RsFunction {
228  public:
238  const int cell,
239  const std::vector<double>& depth,
240  const std::vector<double>& rs)
241  : props_(props)
242  , cell_(cell)
243  , depth_(depth)
244  , rs_(rs)
245  {
246  auto pu = props_.phaseUsage();
247  std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
248  z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
249  z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
250  }
251 
267  double
268  operator()(const double depth,
269  const double press,
270  const double temp,
271  const double sat_gas = 0.0) const
272  {
273  if (sat_gas > 0.0) {
274  return satRs(press, temp);
275  } else {
276  return std::min(satRs(press, temp), linearInterpolation(depth_, rs_, depth));
277  }
278  }
279 
280  private:
281  const BlackoilPropertiesInterface& props_;
282  const int cell_;
283  std::vector<double> depth_;
284  std::vector<double> rs_;
285  double z_[BlackoilPhases::MaxNumPhases];
287 
288  double satRs(const double press, const double temp) const
289  {
290  props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
291  // Rs/Bo is in the gas row and oil column of A_.
292  // 1/Bo is in the oil row and column.
293  // Recall also that it is stored in column-major order.
294  const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
295  const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
296  const int np = props_.numPhases();
297  return A_[np*opos + gpos] / A_[np*opos + opos];
298  }
299  };
300 
301 
307  class RvVD : public RsFunction {
308  public:
318  const int cell,
319  const std::vector<double>& depth,
320  const std::vector<double>& rv)
321  : props_(props)
322  , cell_(cell)
323  , depth_(depth)
324  , rv_(rv)
325  {
326  auto pu = props_.phaseUsage();
327  std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
328  z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
329  z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
330  }
331 
347  double
348  operator()(const double depth,
349  const double press,
350  const double temp,
351  const double sat_oil = 0.0 ) const
352  {
353  if (std::abs(sat_oil) > 1e-16) {
354  return satRv(press, temp);
355  } else {
356  return std::min(satRv(press, temp), linearInterpolation(depth_, rv_, depth));
357  }
358  }
359 
360  private:
361  const BlackoilPropertiesInterface& props_;
362  const int cell_;
363  std::vector<double> depth_;
364  std::vector<double> rv_;
365  double z_[BlackoilPhases::MaxNumPhases];
367 
368  double satRv(const double press, const double temp) const
369  {
370  props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
371  // Rv/Bg is in the oil row and gas column of A_.
372  // 1/Bg is in the gas row and column.
373  // Recall also that it is stored in column-major order.
374  const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
375  const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
376  const int np = props_.numPhases();
377  return A_[np*gpos + opos] / A_[np*gpos + gpos];
378  }
379  };
380 
381 
396  class RsSatAtContact : public RsFunction {
397  public:
406  RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
407  : props_(props), cell_(cell)
408  {
409  auto pu = props_.phaseUsage();
410  std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
411  z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
412  z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
413  rs_sat_contact_ = satRs(p_contact, T_contact);
414  }
415 
431  double
432  operator()(const double /* depth */,
433  const double press,
434  const double temp,
435  const double sat_gas = 0.0) const
436  {
437  if (sat_gas > 0.0) {
438  return satRs(press, temp);
439  } else {
440  return std::min(satRs(press, temp), rs_sat_contact_);
441  }
442  }
443 
444  private:
445  const BlackoilPropertiesInterface& props_;
446  const int cell_;
447  double z_[BlackoilPhases::MaxNumPhases];
448  double rs_sat_contact_;
450 
451  double satRs(const double press, const double temp) const
452  {
453  props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
454  // Rs/Bo is in the gas row and oil column of A_.
455  // 1/Bo is in the oil row and column.
456  // Recall also that it is stored in column-major order.
457  const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
458  const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
459  const int np = props_.numPhases();
460  return A_[np*opos + gpos] / A_[np*opos + opos];
461  }
462  };
463 
464 
479  class RvSatAtContact : public RsFunction {
480  public:
489  RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
490  : props_(props), cell_(cell)
491  {
492  auto pu = props_.phaseUsage();
493  std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
494  z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
495  z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
496  rv_sat_contact_ = satRv(p_contact, T_contact);
497  }
498 
514  double
515  operator()(const double /*depth*/,
516  const double press,
517  const double temp,
518  const double sat_oil = 0.0) const
519  {
520  if (sat_oil > 0.0) {
521  return satRv(press, temp);
522  } else {
523  return std::min(satRv(press, temp), rv_sat_contact_);
524  }
525  }
526 
527  private:
528  const BlackoilPropertiesInterface& props_;
529  const int cell_;
530  double z_[BlackoilPhases::MaxNumPhases];
531  double rv_sat_contact_;
533 
534  double satRv(const double press, const double temp) const
535  {
536  props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
537  // Rv/Bg is in the oil row and gas column of A_.
538  // 1/Bg is in the gas row and column.
539  // Recall also that it is stored in column-major order.
540  const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
541  const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
542  const int np = props_.numPhases();
543  return A_[np*gpos + opos] / A_[np*gpos + gpos];
544  }
545  };
546 
547  } // namespace Miscibility
548 
549 
550 
585  struct EquilRecord {
586  struct {
587  double depth;
588  double press;
589  } main, woc, goc;
592  int N;
593  };
594 
614  template <class DensCalc>
615  class EquilReg {
616  public:
626  EquilReg(const EquilRecord& rec,
627  const DensCalc& density,
628  std::shared_ptr<Miscibility::RsFunction> rs,
629  std::shared_ptr<Miscibility::RsFunction> rv,
630  const PhaseUsage& pu)
631  : rec_ (rec)
632  , density_(density)
633  , rs_ (rs)
634  , rv_ (rv)
635  , pu_ (pu)
636  {
637  }
638 
642  typedef DensCalc CalcDensity;
643 
648 
653 
657  double datum() const { return this->rec_.main.depth; }
658 
662  double pressure() const { return this->rec_.main.press; }
663 
667  double zwoc() const { return this->rec_.woc .depth; }
668 
674  double pcow_woc() const { return this->rec_.woc .press; }
675 
679  double zgoc() const { return this->rec_.goc .depth; }
680 
686  double pcgo_goc() const { return this->rec_.goc .press; }
687 
691  const CalcDensity&
692  densityCalculator() const { return this->density_; }
693 
698  const CalcDissolution&
699  dissolutionCalculator() const { return *this->rs_; }
700 
705  const CalcEvaporation&
706  evaporationCalculator() const { return *this->rv_; }
707 
711  const PhaseUsage&
712  phaseUsage() const { return this->pu_; }
713 
714  private:
715  EquilRecord rec_;
716  DensCalc density_;
717  std::shared_ptr<Miscibility::RsFunction> rs_;
718  std::shared_ptr<Miscibility::RsFunction> rv_;
719  PhaseUsage pu_;
720  };
721 
722 
723 
727  struct PcEq
728  {
730  const int phase,
731  const int cell,
732  const double target_pc)
733  : props_(props),
734  phase_(phase),
735  cell_(cell),
736  target_pc_(target_pc)
737  {
738  std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
739  std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
740  }
741  double operator()(double s) const
742  {
743  s_[phase_] = s;
744  props_.capPress(1, s_, &cell_, pc_, 0);
745  return pc_[phase_] - target_pc_;
746  }
747  private:
748  const BlackoilPropertiesInterface& props_;
749  const int phase_;
750  const int cell_;
751  const int target_pc_;
752  mutable double s_[BlackoilPhases::MaxNumPhases];
753  mutable double pc_[BlackoilPhases::MaxNumPhases];
754  };
755 
756 
757 
760  inline double satFromPc(const BlackoilPropertiesInterface& props,
761  const int phase,
762  const int cell,
763  const double target_pc,
764  const bool increasing = false)
765  {
766  // Find minimum and maximum saturations.
767  double sminarr[BlackoilPhases::MaxNumPhases];
768  double smaxarr[BlackoilPhases::MaxNumPhases];
769  props.satRange(1, &cell, sminarr, smaxarr);
770  const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
771  const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
772 
773  // Create the equation f(s) = pc(s) - target_pc
774  const PcEq f(props, phase, cell, target_pc);
775  const double f0 = f(s0);
776  const double f1 = f(s1);
777 
778  if (f0 <= 0.0) {
779  return s0;
780  } else if (f1 > 0.0) {
781  return s1;
782  } else {
783  const int max_iter = 30;
784  const double tol = 1e-6;
785  int iter_used = -1;
786  typedef RegulaFalsi<ThrowOnError> ScalarSolver;
787  const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
788  return sol;
789  }
790  }
791 
792 
796  struct PcEqSum
797  {
799  const int phase1,
800  const int phase2,
801  const int cell,
802  const double target_pc)
803  : props_(props),
804  phase1_(phase1),
805  phase2_(phase2),
806  cell_(cell),
807  target_pc_(target_pc)
808  {
809  std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
810  std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
811  }
812  double operator()(double s) const
813  {
814  s_[phase1_] = s;
815  s_[phase2_] = 1.0 - s;
816  props_.capPress(1, s_, &cell_, pc_, 0);
817  return pc_[phase1_] + pc_[phase2_] - target_pc_;
818  }
819  private:
820  const BlackoilPropertiesInterface& props_;
821  const int phase1_;
822  const int phase2_;
823  const int cell_;
824  const int target_pc_;
825  mutable double s_[BlackoilPhases::MaxNumPhases];
826  mutable double pc_[BlackoilPhases::MaxNumPhases];
827  };
828 
829 
830 
831 
835  inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
836  const int phase1,
837  const int phase2,
838  const int cell,
839  const double target_pc)
840  {
841  // Find minimum and maximum saturations.
842  double sminarr[BlackoilPhases::MaxNumPhases];
843  double smaxarr[BlackoilPhases::MaxNumPhases];
844  props.satRange(1, &cell, sminarr, smaxarr);
845  const double smin = sminarr[phase1];
846  const double smax = smaxarr[phase1];
847 
848  // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
849  const PcEqSum f(props, phase1, phase2, cell, target_pc);
850  const double f0 = f(smin);
851  const double f1 = f(smax);
852  if (f0 <= 0.0) {
853  return smin;
854  } else if (f1 > 0.0) {
855  return smax;
856  } else {
857  const int max_iter = 30;
858  const double tol = 1e-6;
859  int iter_used = -1;
860  typedef RegulaFalsi<ThrowOnError> ScalarSolver;
861  const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
862  return sol;
863  }
864  }
865 
867  inline double satFromDepth(const BlackoilPropertiesInterface& props,
868  const double cellDepth,
869  const double contactDepth,
870  const int phase,
871  const int cell,
872  const bool increasing = false)
873  {
874  // Find minimum and maximum saturations.
875  double sminarr[BlackoilPhases::MaxNumPhases];
876  double smaxarr[BlackoilPhases::MaxNumPhases];
877  props.satRange(1, &cell, sminarr, smaxarr);
878  const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
879  const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
880 
881  if (cellDepth < contactDepth){
882  return s0;
883  } else {
884  return s1;
885  }
886 
887  }
888 
890  inline bool isConstPc(const BlackoilPropertiesInterface& props,
891  const int phase,
892  const int cell)
893  {
894  // Find minimum and maximum saturations.
895  double sminarr[BlackoilPhases::MaxNumPhases];
896  double smaxarr[BlackoilPhases::MaxNumPhases];
897  props.satRange(1, &cell, sminarr, smaxarr);
898 
899  // Create the equation f(s) = pc(s);
900  const PcEq f(props, phase, cell, 0);
901  const double f0 = f(sminarr[phase]);
902  const double f1 = f(smaxarr[phase]);
903  return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
904  }
905 
906  } // namespace Equil
907 } // namespace Opm
908 
909 
910 #endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
EquilReg(const EquilRecord &rec, const DensCalc &density, std::shared_ptr< Miscibility::RsFunction > rs, std::shared_ptr< Miscibility::RsFunction > rv, const PhaseUsage &pu)
Definition: EquilibrationHelpers.hpp:626
Definition: EquilibrationHelpers.hpp:396
double linearInterpolation(const std::vector< double > &xv, const std::vector< double > &yv, double x)
Definition: linearInterpolation.hpp:64
virtual double operator()(const double depth, const double press, const double temp, const double sat=0.0) const =0
Definition: EquilibrationHelpers.hpp:479
int live_oil_table_index
Definition: EquilibrationHelpers.hpp:590
struct Opm::Equil::EquilRecord::@4 woc
Definition: AnisotropicEikonal.hpp:43
Definition: BlackoilPhases.hpp:32
DensityCalculator(const BlackoilPropertiesInterface &props, const int c)
Definition: EquilibrationHelpers.hpp:112
double operator()(const double, const double, const double, const double=0.0) const
Definition: EquilibrationHelpers.hpp:212
virtual void matrix(const int n, const double *p, const double *T, const double *z, const int *cells, double *A, double *dAdp) const =0
double operator()(const double, const double press, const double temp, const double sat_oil=0.0) const
Definition: EquilibrationHelpers.hpp:515
struct Opm::Equil::EquilRecord::@4 main
double zgoc() const
Definition: EquilibrationHelpers.hpp:679
const CalcDissolution & dissolutionCalculator() const
Definition: EquilibrationHelpers.hpp:699
virtual void capPress(const int n, const double *s, const int *cells, double *pc, double *dpcds) const =0
double press
Definition: EquilibrationHelpers.hpp:588
Definition: RootFinders.hpp:103
double satFromDepth(const BlackoilPropertiesInterface &props, const double cellDepth, const double contactDepth, const int phase, const int cell, const bool increasing=false)
Compute saturation from depth. Used for constant capillary pressure function.
Definition: EquilibrationHelpers.hpp:867
bool isConstPc(const BlackoilPropertiesInterface &props, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers.hpp:890
virtual int numPhases() const =0
const CalcEvaporation & evaporationCalculator() const
Definition: EquilibrationHelpers.hpp:706
double operator()(const double depth, const double press, const double temp, const double sat_oil=0.0) const
Definition: EquilibrationHelpers.hpp:348
DensCalc CalcDensity
Definition: EquilibrationHelpers.hpp:642
Definition: BlackoilPropertiesInterface.hpp:37
int phase_pos[MaxNumPhases]
Definition: BlackoilPhases.hpp:40
Definition: EquilibrationHelpers.hpp:91
double operator()(const double depth, const double press, const double temp, const double sat_gas=0.0) const
Definition: EquilibrationHelpers.hpp:268
virtual void satRange(const int n, const int *cells, double *smin, double *smax) const =0
double zwoc() const
Definition: EquilibrationHelpers.hpp:667
double datum() const
Definition: EquilibrationHelpers.hpp:657
virtual PhaseUsage phaseUsage() const =0
double satFromPc(const BlackoilPropertiesInterface &props, const int phase, const int cell, const double target_pc, const bool increasing=false)
Definition: EquilibrationHelpers.hpp:760
Miscibility::RsFunction CalcDissolution
Definition: EquilibrationHelpers.hpp:647
PcEq(const BlackoilPropertiesInterface &props, const int phase, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:729
double operator()(double s) const
Definition: EquilibrationHelpers.hpp:812
int wet_gas_table_index
Definition: EquilibrationHelpers.hpp:591
double pcgo_goc() const
Definition: EquilibrationHelpers.hpp:686
Definition: EquilibrationHelpers.hpp:796
struct Opm::Equil::EquilRecord::@4 goc
Definition: EquilibrationHelpers.hpp:193
double pcow_woc() const
Definition: EquilibrationHelpers.hpp:674
int N
Definition: EquilibrationHelpers.hpp:592
double operator()(double s) const
Definition: EquilibrationHelpers.hpp:741
Definition: EquilibrationHelpers.hpp:585
RsSatAtContact(const BlackoilPropertiesInterface &props, const int cell, const double p_contact, const double T_contact)
Definition: EquilibrationHelpers.hpp:406
const PhaseUsage & phaseUsage() const
Definition: EquilibrationHelpers.hpp:712
Definition: EquilibrationHelpers.hpp:615
RvVD(const BlackoilPropertiesInterface &props, const int cell, const std::vector< double > &depth, const std::vector< double > &rv)
Definition: EquilibrationHelpers.hpp:317
Definition: BlackoilPhases.hpp:32
double pressure() const
Definition: EquilibrationHelpers.hpp:662
double operator()(const double, const double press, const double temp, const double sat_gas=0.0) const
Definition: EquilibrationHelpers.hpp:432
static const int MaxNumPhases
Definition: BlackoilPhases.hpp:30
RvSatAtContact(const BlackoilPropertiesInterface &props, const int cell, const double p_contact, const double T_contact)
Definition: EquilibrationHelpers.hpp:489
RsVD(const BlackoilPropertiesInterface &props, const int cell, const std::vector< double > &depth, const std::vector< double > &rs)
Definition: EquilibrationHelpers.hpp:237
const CalcDensity & densityCalculator() const
Definition: EquilibrationHelpers.hpp:692
PcEqSum(const BlackoilPropertiesInterface &props, const int phase1, const int phase2, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:798
Definition: EquilibrationHelpers.hpp:307
double depth
Definition: EquilibrationHelpers.hpp:587
Definition: BlackoilPhases.hpp:36
Definition: EquilibrationHelpers.hpp:227
Miscibility::RsFunction CalcEvaporation
Definition: EquilibrationHelpers.hpp:652
Definition: EquilibrationHelpers.hpp:165
Definition: EquilibrationHelpers.hpp:727
double satFromSumOfPcs(const BlackoilPropertiesInterface &props, const int phase1, const int phase2, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:835
std::vector< double > operator()(const double p, const double T, const std::vector< double > &z) const
Definition: EquilibrationHelpers.hpp:132