initStateEquil_impl.hpp
1 /*
2  Copyright 2014 SINTEF ICT, Applied Mathematics.
3  Copyright 2015 Dr. Blatt - HPC-Simulation-Software & Services
4  Copyright 2015 NTNU
5 
6  This file is part of the Open Porous Media project (OPM).
7 
8  OPM is free software: you can redistribute it and/or modify
9  it under the terms of the GNU General Public License as published by
10  the Free Software Foundation, either version 3 of the License, or
11  (at your option) any later version.
12 
13  OPM is distributed in the hope that it will be useful,
14  but WITHOUT ANY WARRANTY; without even the implied warranty of
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16  GNU General Public License for more details.
17 
18  You should have received a copy of the GNU General Public License
19  along with OPM. If not, see <http://www.gnu.org/licenses/>.
20 */
21 
22 #ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
23 #define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
24 
25 #include <opm/core/grid.h>
26 #include <opm/core/grid/GridHelpers.hpp>
27 #include <opm/core/props/BlackoilPhases.hpp>
29 
30 #include <cassert>
31 #include <cmath>
32 #include <functional>
33 #include <vector>
34 
35 namespace Opm
36 {
37  namespace Details {
38  template <class RHS>
39  class RK4IVP {
40  public:
41  RK4IVP(const RHS& f ,
42  const std::array<double,2>& span,
43  const double y0 ,
44  const int N )
45  : N_(N)
46  , span_(span)
47  {
48  const double h = stepsize();
49  const double h2 = h / 2;
50  const double h6 = h / 6;
51 
52  y_.reserve(N + 1);
53  f_.reserve(N + 1);
54 
55  y_.push_back(y0);
56  f_.push_back(f(span_[0], y0));
57 
58  for (int i = 0; i < N; ++i) {
59  const double x = span_[0] + i*h;
60  const double y = y_.back();
61 
62  const double k1 = f_[i];
63  const double k2 = f(x + h2, y + h2*k1);
64  const double k3 = f(x + h2, y + h2*k2);
65  const double k4 = f(x + h , y + h*k3);
66 
67  y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
68  f_.push_back(f(x + h, y_.back()));
69  }
70 
71  assert (y_.size() == std::vector<double>::size_type(N + 1));
72  }
73 
74  double
75  operator()(const double x) const
76  {
77  // Dense output (O(h**3)) according to Shampine
78  // (Hermite interpolation)
79  const double h = stepsize();
80  int i = (x - span_[0]) / h;
81  const double t = (x - (span_[0] + i*h)) / h;
82 
83  // Crude handling of evaluation point outside "span_";
84  if (i < 0) { i = 0; }
85  if (N_ <= i) { i = N_ - 1; }
86 
87  const double y0 = y_[i], y1 = y_[i + 1];
88  const double f0 = f_[i], f1 = f_[i + 1];
89 
90  double u = (1 - 2*t) * (y1 - y0);
91  u += h * ((t - 1)*f0 + t*f1);
92  u *= t * (t - 1);
93  u += (1 - t)*y0 + t*y1;
94 
95  return u;
96  }
97 
98  private:
99  int N_;
100  std::array<double,2> span_;
101  std::vector<double> y_;
102  std::vector<double> f_;
103 
104  double
105  stepsize() const { return (span_[1] - span_[0]) / N_; }
106  };
107 
108  namespace PhasePressODE {
109  template <class Density>
110  class Water {
111  public:
112  Water(const double temp,
113  const Density& rho,
114  const int np,
115  const int ix,
116  const double norm_grav)
117  : temp_(temp)
118  , rho_(rho)
119  , svol_(np, 0)
120  , ix_(ix)
121  , g_(norm_grav)
122  {
123  svol_[ix_] = 1.0;
124  }
125 
126  double
127  operator()(const double /* depth */,
128  const double press) const
129  {
130  return this->density(press) * g_;
131  }
132 
133  private:
134  const double temp_;
135  const Density& rho_;
136  std::vector<double> svol_;
137  const int ix_;
138  const double g_;
139 
140  double
141  density(const double press) const
142  {
143  const std::vector<double>& rho = rho_(press, temp_, svol_);
144 
145  return rho[ix_];
146  }
147  };
148 
149  template <class Density, class RS>
150  class Oil {
151  public:
152  Oil(const double temp,
153  const Density& rho,
154  const RS& rs,
155  const int np,
156  const int oix,
157  const int gix,
158  const double norm_grav)
159  : temp_(temp)
160  , rho_(rho)
161  , rs_(rs)
162  , svol_(np, 0)
163  , oix_(oix)
164  , gix_(gix)
165  , g_(norm_grav)
166  {
167  svol_[oix_] = 1.0;
168  }
169 
170  double
171  operator()(const double depth,
172  const double press) const
173  {
174  return this->density(depth, press) * g_;
175  }
176 
177  private:
178  const double temp_;
179  const Density& rho_;
180  const RS& rs_;
181  mutable std::vector<double> svol_;
182  const int oix_;
183  const int gix_;
184  const double g_;
185 
186  double
187  density(const double depth,
188  const double press) const
189  {
190  if (gix_ >= 0) {
191  svol_[gix_] = rs_(depth, press, temp_);
192  }
193 
194  const std::vector<double>& rho = rho_(press, temp_, svol_);
195  return rho[oix_];
196  }
197  };
198 
199  template <class Density, class RV>
200  class Gas {
201  public:
202  Gas(const double temp,
203  const Density& rho,
204  const RV& rv,
205  const int np,
206  const int gix,
207  const int oix,
208  const double norm_grav)
209  : temp_(temp)
210  , rho_(rho)
211  , rv_(rv)
212  , svol_(np, 0)
213  , gix_(gix)
214  , oix_(oix)
215  , g_(norm_grav)
216  {
217  svol_[gix_] = 1.0;
218  }
219 
220  double
221  operator()(const double depth,
222  const double press) const
223  {
224  return this->density(depth, press) * g_;
225  }
226 
227  private:
228  const double temp_;
229  const Density& rho_;
230  const RV& rv_;
231  mutable std::vector<double> svol_;
232  const int gix_;
233  const int oix_;
234  const double g_;
235 
236  double
237  density(const double depth,
238  const double press) const
239  {
240  if (oix_ >= 0) {
241  svol_[oix_] = rv_(depth, press, temp_);
242  }
243 
244  const std::vector<double>& rho = rho_(press, temp_, svol_);
245  return rho[gix_];
246  }
247  };
248  } // namespace PhasePressODE
249 
250  namespace PhaseUsed {
251  inline bool
252  water(const PhaseUsage& pu)
253  {
254  return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
255  }
256 
257  inline bool
258  oil(const PhaseUsage& pu)
259  {
260  return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
261  }
262 
263  inline bool
264  gas(const PhaseUsage& pu)
265  {
266  return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
267  }
268  } // namespace PhaseUsed
269 
270  namespace PhaseIndex {
271  inline int
272  water(const PhaseUsage& pu)
273  {
274  int i = -1;
275  if (PhaseUsed::water(pu)) {
276  i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
277  }
278 
279  return i;
280  }
281 
282  inline int
283  oil(const PhaseUsage& pu)
284  {
285  int i = -1;
286  if (PhaseUsed::oil(pu)) {
287  i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
288  }
289 
290  return i;
291  }
292 
293  inline int
294  gas(const PhaseUsage& pu)
295  {
296  int i = -1;
297  if (PhaseUsed::gas(pu)) {
298  i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
299  }
300 
301  return i;
302  }
303  } // namespace PhaseIndex
304 
305  namespace PhasePressure {
306  template <class Grid,
307  class PressFunction,
308  class CellRange>
309  void
310  assign(const Grid& G ,
311  const std::array<PressFunction, 2>& f ,
312  const double split,
313  const CellRange& cells,
314  std::vector<double>& p )
315  {
316 
317  enum { up = 0, down = 1 };
318 
319  std::vector<double>::size_type c = 0;
320  for (typename CellRange::const_iterator
321  ci = cells.begin(), ce = cells.end();
322  ci != ce; ++ci, ++c)
323  {
324  assert (c < p.size());
325 
326  const double z = UgGridHelpers::cellCenterDepth(G, *ci);
327  p[c] = (z < split) ? f[up](z) : f[down](z);
328  }
329  }
330 
331  template <class Grid,
332  class Region,
333  class CellRange>
334  void
335  water(const Grid& G ,
336  const Region& reg ,
337  const std::array<double,2>& span ,
338  const double grav ,
339  double& po_woc,
340  const CellRange& cells ,
341  std::vector<double>& press )
342  {
343  using PhasePressODE::Water;
344  typedef Water<typename Region::CalcDensity> ODE;
345 
346  const PhaseUsage& pu = reg.phaseUsage();
347 
348  const int wix = PhaseIndex::water(pu);
349  const double T = 273.15 + 20; // standard temperature for now
350  ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
351 
352  double z0;
353  double p0;
354  if (reg.datum() > reg.zwoc()) {//Datum in water zone
355  z0 = reg.datum();
356  p0 = reg.pressure();
357  } else {
358  z0 = reg.zwoc();
359  p0 = po_woc - reg.pcow_woc(); // Water pressure at contact
360  }
361 
362  std::array<double,2> up = {{ z0, span[0] }};
363  std::array<double,2> down = {{ z0, span[1] }};
364 
365  typedef Details::RK4IVP<ODE> WPress;
366  std::array<WPress,2> wpress = {
367  {
368  WPress(drho, up , p0, 2000)
369  ,
370  WPress(drho, down, p0, 2000)
371  }
372  };
373 
374  assign(G, wpress, z0, cells, press);
375 
376  if (reg.datum() > reg.zwoc()) {
377  // Return oil pressure at contact
378  po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc();
379  }
380  }
381 
382  template <class Grid,
383  class Region,
384  class CellRange>
385  void
386  oil(const Grid& G ,
387  const Region& reg ,
388  const std::array<double,2>& span ,
389  const double grav ,
390  const CellRange& cells ,
391  std::vector<double>& press ,
392  double& po_woc,
393  double& po_goc)
394  {
395  using PhasePressODE::Oil;
396  typedef Oil<typename Region::CalcDensity,
397  typename Region::CalcDissolution> ODE;
398 
399  const PhaseUsage& pu = reg.phaseUsage();
400 
401  const int oix = PhaseIndex::oil(pu);
402  const int gix = PhaseIndex::gas(pu);
403  const double T = 273.15 + 20; // standard temperature for now
404  ODE drho(T,
405  reg.densityCalculator(),
406  reg.dissolutionCalculator(),
407  pu.num_phases, oix, gix, grav);
408 
409  double z0;
410  double p0;
411  if (reg.datum() > reg.zwoc()) {//Datum in water zone, po_woc given
412  z0 = reg.zwoc();
413  p0 = po_woc;
414  } else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, po_goc given
415  z0 = reg.zgoc();
416  p0 = po_goc;
417  } else { //Datum in oil zone
418  z0 = reg.datum();
419  p0 = reg.pressure();
420  }
421 
422  std::array<double,2> up = {{ z0, span[0] }};
423  std::array<double,2> down = {{ z0, span[1] }};
424 
425  typedef Details::RK4IVP<ODE> OPress;
426  std::array<OPress,2> opress = {
427  {
428  OPress(drho, up , p0, 2000)
429  ,
430  OPress(drho, down, p0, 2000)
431  }
432  };
433 
434  assign(G, opress, z0, cells, press);
435 
436  const double woc = reg.zwoc();
437  if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum
438  else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum
439  else { po_woc = p0; } // WOC *at* datum
440 
441  const double goc = reg.zgoc();
442  if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum
443  else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum
444  else { po_goc = p0; } // GOC *at* datum
445  }
446 
447  template <class Grid,
448  class Region,
449  class CellRange>
450  void
451  gas(const Grid& G ,
452  const Region& reg ,
453  const std::array<double,2>& span ,
454  const double grav ,
455  double& po_goc,
456  const CellRange& cells ,
457  std::vector<double>& press )
458  {
459  using PhasePressODE::Gas;
460  typedef Gas<typename Region::CalcDensity,
461  typename Region::CalcEvaporation> ODE;
462 
463  const PhaseUsage& pu = reg.phaseUsage();
464 
465  const int gix = PhaseIndex::gas(pu);
466  const int oix = PhaseIndex::oil(pu);
467 
468  const double T = 273.15 + 20; // standard temperature for now
469  ODE drho(T,
470  reg.densityCalculator(),
471  reg.evaporationCalculator(),
472  pu.num_phases, gix, oix, grav);
473 
474  double z0;
475  double p0;
476  if (reg.datum() < reg.zgoc()) {//Datum in gas zone
477  z0 = reg.datum();
478  p0 = reg.pressure();
479  } else {
480  z0 = reg.zgoc();
481  p0 = po_goc + reg.pcgo_goc(); // Gas pressure at contact
482  }
483 
484  std::array<double,2> up = {{ z0, span[0] }};
485  std::array<double,2> down = {{ z0, span[1] }};
486 
487  typedef Details::RK4IVP<ODE> GPress;
488  std::array<GPress,2> gpress = {
489  {
490  GPress(drho, up , p0, 2000)
491  ,
492  GPress(drho, down, p0, 2000)
493  }
494  };
495 
496  assign(G, gpress, z0, cells, press);
497 
498  if (reg.datum() < reg.zgoc()) {
499  // Return oil pressure at contact
500  po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc();
501  }
502  }
503  } // namespace PhasePressure
504 
505  template <class Grid,
506  class Region,
507  class CellRange>
508  void
509  equilibrateOWG(const Grid& G,
510  const Region& reg,
511  const double grav,
512  const std::array<double,2>& span,
513  const CellRange& cells,
514  std::vector< std::vector<double> >& press)
515  {
516  const PhaseUsage& pu = reg.phaseUsage();
517 
518  if (reg.datum() > reg.zwoc()) { // Datum in water zone
519  double po_woc = -1;
520  double po_goc = -1;
521 
522  if (PhaseUsed::water(pu)) {
523  const int wix = PhaseIndex::water(pu);
524  PhasePressure::water(G, reg, span, grav, po_woc,
525  cells, press[ wix ]);
526  }
527 
528  if (PhaseUsed::oil(pu)) {
529  const int oix = PhaseIndex::oil(pu);
530  PhasePressure::oil(G, reg, span, grav, cells,
531  press[ oix ], po_woc, po_goc);
532  }
533 
534  if (PhaseUsed::gas(pu)) {
535  const int gix = PhaseIndex::gas(pu);
536  PhasePressure::gas(G, reg, span, grav, po_goc,
537  cells, press[ gix ]);
538  }
539  } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
540  double po_woc = -1;
541  double po_goc = -1;
542 
543  if (PhaseUsed::gas(pu)) {
544  const int gix = PhaseIndex::gas(pu);
545  PhasePressure::gas(G, reg, span, grav, po_goc,
546  cells, press[ gix ]);
547  }
548 
549  if (PhaseUsed::oil(pu)) {
550  const int oix = PhaseIndex::oil(pu);
551  PhasePressure::oil(G, reg, span, grav, cells,
552  press[ oix ], po_woc, po_goc);
553  }
554 
555  if (PhaseUsed::water(pu)) {
556  const int wix = PhaseIndex::water(pu);
557  PhasePressure::water(G, reg, span, grav, po_woc,
558  cells, press[ wix ]);
559  }
560  } else { // Datum in oil zone
561  double po_woc = -1;
562  double po_goc = -1;
563 
564  if (PhaseUsed::oil(pu)) {
565  const int oix = PhaseIndex::oil(pu);
566  PhasePressure::oil(G, reg, span, grav, cells,
567  press[ oix ], po_woc, po_goc);
568  }
569 
570  if (PhaseUsed::water(pu)) {
571  const int wix = PhaseIndex::water(pu);
572  PhasePressure::water(G, reg, span, grav, po_woc,
573  cells, press[ wix ]);
574  }
575 
576  if (PhaseUsed::gas(pu)) {
577  const int gix = PhaseIndex::gas(pu);
578  PhasePressure::gas(G, reg, span, grav, po_goc,
579  cells, press[ gix ]);
580  }
581  }
582  }
583  } // namespace Details
584 
585 
586  namespace EQUIL {
587 
588 
589  template <class Grid,
590  class Region,
591  class CellRange>
592  std::vector< std::vector<double> >
593  phasePressures(const Grid& G,
594  const Region& reg,
595  const CellRange& cells,
596  const double grav)
597  {
598  std::array<double,2> span =
599  {{ std::numeric_limits<double>::max() ,
600  -std::numeric_limits<double>::max() }}; // Symm. about 0.
601 
602  int ncell = 0;
603  {
604  // This code is only supported in three space dimensions
605  assert (UgGridHelpers::dimensions(G) == 3);
606 
607  const int nd = UgGridHelpers::dimensions(G);
608 
609  // Define vertical span as
610  //
611  // [minimum(node depth(cells)), maximum(node depth(cells))]
612  //
613  // Note: We use a sledgehammer approach--looping all
614  // the nodes of all the faces of all the 'cells'--to
615  // compute those bounds. This necessarily entails
616  // visiting some nodes (and faces) multiple times.
617  //
618  // Note: The implementation of 'RK4IVP<>' implicitly
619  // imposes the requirement that cell centroids are all
620  // within this vertical span. That requirement is not
621  // checked.
622  auto cell2Faces = UgGridHelpers::cell2Faces(G);
623  auto faceVertices = UgGridHelpers::face2Vertices(G);
624 
625  for (typename CellRange::const_iterator
626  ci = cells.begin(), ce = cells.end();
627  ci != ce; ++ci, ++ncell)
628  {
629  for (auto fi=cell2Faces[*ci].begin(),
630  fe=cell2Faces[*ci].end();
631  fi != fe;
632  ++fi)
633  {
634  for (auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end();
635  i != e; ++i)
636  {
637  const double z = UgGridHelpers::vertexCoordinates(G, *i)[nd-1];
638 
639  if (z < span[0]) { span[0] = z; }
640  if (z > span[1]) { span[1] = z; }
641  }
642  }
643  }
644  }
645  const int np = reg.phaseUsage().num_phases;
646 
647  typedef std::vector<double> pval;
648  std::vector<pval> press(np, pval(ncell, 0.0));
649 
650  const double zwoc = reg.zwoc ();
651  const double zgoc = reg.zgoc ();
652 
653  // make sure goc and woc is within the span for the phase pressure calculation
654  span[0] = std::min(span[0],zgoc);
655  span[1] = std::max(span[1],zwoc);
656 
657  Details::equilibrateOWG(G, reg, grav, span, cells, press);
658 
659  return press;
660  }
661 
662  template <class Grid,
663  class Region,
664  class CellRange>
665  std::vector<double>
666  temperature(const Grid& /* G */,
667  const Region& /* reg */,
668  const CellRange& cells)
669  {
670  // use the standard temperature for everything for now
671  return std::vector<double>(cells.size(), 273.15 + 20.0);
672  }
673 
674  template <class Grid, class Region, class CellRange>
675  std::vector< std::vector<double> >
676  phaseSaturations(const Grid& G,
677  const Region& reg,
678  const CellRange& cells,
679  BlackoilPropertiesInterface& props,
680  const std::vector<double> swat_init,
681  std::vector< std::vector<double> >& phase_pressures)
682  {
683  if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
684  OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
685  }
686 
687  std::vector< std::vector<double> > phase_saturations = phase_pressures; // Just to get the right size.
688  double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
689  double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
690 
691  const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
692  const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
693  const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
694  const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
695  const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
696  std::vector<double>::size_type local_index = 0;
697  for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
698  const int cell = *ci;
699  props.satRange(1, &cell, smin, smax);
700  // Find saturations from pressure differences by
701  // inverting capillary pressure functions.
702  double sw = 0.0;
703  if (water) {
704  if (isConstPc(props,waterpos,cell)){
705  const double cellDepth = UgGridHelpers::cellCenterDepth(G,
706  cell);
707  sw = satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,false);
708  phase_saturations[waterpos][local_index] = sw;
709  }
710  else{
711  const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
712  if (swat_init.empty()) { // Invert Pc to find sw
713  sw = satFromPc(props, waterpos, cell, pcov);
714  phase_saturations[waterpos][local_index] = sw;
715  } else { // Scale Pc to reflect imposed sw
716  sw = swat_init[cell];
717  props.swatInitScaling(cell, pcov, sw);
718  phase_saturations[waterpos][local_index] = sw;
719  }
720  }
721  }
722  double sg = 0.0;
723  if (gas) {
724  if (isConstPc(props,gaspos,cell)){
725  const double cellDepth = UgGridHelpers::cellCenterDepth(G,
726  cell);
727  sg = satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,true);
728  phase_saturations[gaspos][local_index] = sg;
729  }
730  else{
731  // Note that pcog is defined to be (pg - po), not (po - pg).
732  const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
733  const double increasing = true; // pcog(sg) expected to be increasing function
734  sg = satFromPc(props, gaspos, cell, pcog, increasing);
735  phase_saturations[gaspos][local_index] = sg;
736  }
737  }
738  if (gas && water && (sg + sw > 1.0)) {
739  // Overlapping gas-oil and oil-water transition
740  // zones can lead to unphysical saturations when
741  // treated as above. Must recalculate using gas-water
742  // capillary pressure.
743  const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
744  if (! swat_init.empty()) {
745  // Re-scale Pc to reflect imposed sw for vanishing oil phase.
746  // This seems consistent with ecl, and fails to honour
747  // swat_init in case of non-trivial gas-oil cap pressure.
748  props.swatInitScaling(cell, pcgw, sw);
749  }
750  sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
751  sg = 1.0 - sw;
752  phase_saturations[waterpos][local_index] = sw;
753  phase_saturations[gaspos][local_index] = sg;
754  // Adjust oil pressure according to gas saturation and cap pressure
755  double pc[BlackoilPhases::MaxNumPhases];
756  double sat[BlackoilPhases::MaxNumPhases];
757  sat[waterpos] = sw;
758  sat[gaspos] = sg;
759  sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
760  props.capPress(1, sat, &cell, pc, 0);
761  phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
762  }
763  phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
764 
765  // Adjust phase pressures for max and min saturation ...
766  double pc[BlackoilPhases::MaxNumPhases];
767  double sat[BlackoilPhases::MaxNumPhases];
768  double threshold_sat = 1.0e-6;
769 
770  sat[oilpos] = 1.0;
771  if (water) {
772  sat[waterpos] = smax[waterpos];
773  sat[oilpos] -= sat[waterpos];
774  }
775  if (gas) {
776  sat[gaspos] = smax[gaspos];
777  sat[oilpos] -= sat[gaspos];
778  }
779  if (water && sw > smax[waterpos]-threshold_sat ) {
780  sat[waterpos] = smax[waterpos];
781  props.capPress(1, sat, &cell, pc, 0);
782  phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
783  } else if (gas && sg > smax[gaspos]-threshold_sat) {
784  sat[gaspos] = smax[gaspos];
785  props.capPress(1, sat, &cell, pc, 0);
786  phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
787  }
788  if (gas && sg < smin[gaspos]+threshold_sat) {
789  sat[gaspos] = smin[gaspos];
790  props.capPress(1, sat, &cell, pc, 0);
791  phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
792  }
793  if (water && sw < smin[waterpos]+threshold_sat) {
794  sat[waterpos] = smin[waterpos];
795  props.capPress(1, sat, &cell, pc, 0);
796  phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
797  }
798  }
799  return phase_saturations;
800  }
801 
802 
821  template <class Grid, class CellRangeType>
822  std::vector<double> computeRs(const Grid& grid,
823  const CellRangeType& cells,
824  const std::vector<double> oil_pressure,
825  const std::vector<double>& temperature,
826  const Miscibility::RsFunction& rs_func,
827  const std::vector<double> gas_saturation)
828  {
829  assert(UgGridHelpers::dimensions(grid) == 3);
830  std::vector<double> rs(cells.size());
831  int count = 0;
832  for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
833  const double depth = UgGridHelpers::cellCenterDepth(grid, *it);
834  rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
835  }
836  return rs;
837  }
838 
839  } // namespace Equil
840 
841 
842  namespace Details
843  {
847  inline std::vector<double>
848  convertSats(const std::vector< std::vector<double> >& sat)
849  {
850  const auto np = sat.size();
851  const auto nc = sat[0].size();
852 
853  std::vector<double> s(np * nc);
854 
855  for (decltype(sat.size()) p = 0; p < np; ++p) {
856  const auto& sat_p = sat[p];
857  double* sp = & s[0*nc + p];
858 
859  for (decltype(sat[0].size()) c = 0;
860  c < nc; ++c, sp += np)
861  {
862  *sp = sat_p[c];
863  }
864  }
865 
866  return s;
867  }
868  } // namespace Details
869 
870 
888  template<class Grid>
889  void initStateEquil(const Grid& grid,
890  BlackoilPropertiesFromDeck& props,
891  const Opm::Deck& deck,
892  const Opm::EclipseState& eclipseState,
893  const double gravity,
894  BlackoilState& state,
895  bool applySwatinit = true)
896  {
897 
898  typedef EQUIL::DeckDependent::InitialStateComputer ISC;
899  //Check for presence of kw SWATINIT
900  std::vector<double> swat_init = {};
901  if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatinit) {
902  const std::vector<double>& swat_init_ecl = eclipseState.
903  get3DProperties().getDoubleGridProperty("SWATINIT").getData();
904  const int nc = UgGridHelpers::numCells(grid);
905  swat_init.resize(nc);
906  const int* gc = UgGridHelpers::globalCell(grid);
907  for (int c = 0; c < nc; ++c) {
908  const int deck_pos = (gc == NULL) ? c : gc[c];
909  swat_init[c] = swat_init_ecl[deck_pos];
910  }
911  }
912 
913  ISC isc(props, deck, eclipseState, grid, gravity, swat_init);
914  const auto pu = props.phaseUsage();
915  const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
916  ? pu.phase_pos[BlackoilPhases::Liquid]
917  : pu.phase_pos[BlackoilPhases::Aqua];
918  state.pressure() = isc.press()[ref_phase];
919  state.saturation() = Details::convertSats(isc.saturation());
920  state.gasoilratio() = isc.rs();
921  state.rv() = isc.rv();
922 
923  initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
924  }
925 
926 
927 
928 } // namespace Opm
929 
930 #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
void initStateEquil(const Grid &grid, const BlackoilPropertiesInterface &props, const Opm::Deck &deck, const Opm::EclipseState &eclipseState, const double gravity, BlackoilState &state, bool applySwatInit=true)
Compute initial state by an equilibration procedure.
double satFromPc(const BlackoilPropertiesInterface &props, const int phase, const int cell, const double target_pc, const bool increasing=false)
Compute saturation of some phase corresponding to a given capillary pressure.
Definition: EquilibrationHelpers.hpp:714
Definition: initStateEquil_impl.hpp:110
Definition: initStateEquil_impl.hpp:39
Definition: AnisotropicEikonal.cpp:446
Definition: initStateEquil_impl.hpp:150
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:821
Functions for initializing a reservoir state.
Definition: BlackoilPhases.hpp:36
bool isConstPc(const BlackoilPropertiesInterface &props, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers.hpp:844
Definition: initStateEquil_impl.hpp:200
Base class for phase mixing functions.
Definition: EquilibrationHelpers.hpp:165
double satFromSumOfPcs(const BlackoilPropertiesInterface &props, const int phase1, const int phase2, const int cell, const double target_pc)
Compute saturation of some phase corresponding to a given capillary pressure, where the capillary pre...
Definition: EquilibrationHelpers.hpp:789
std::vector< double > computeRs(const Grid &grid, const CellRangeType &cells, const std::vector< double > oil_pressure, const std::vector< double > &temperature, const Miscibility::RsFunction &rs_func, const std::vector< double > gas_saturation)
Compute initial Rs values.
Definition: initStateEquil_impl.hpp:822