initStateEquil_impl.hpp
Go to the documentation of this file.
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>
28 #include <opm/core/simulator/initState.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  const int nd = UgGridHelpers::dimensions(G);
317 
318  enum { up = 0, down = 1 };
319 
320  std::vector<double>::size_type c = 0;
321  for (typename CellRange::const_iterator
322  ci = cells.begin(), ce = cells.end();
323  ci != ce; ++ci, ++c)
324  {
325  assert (c < p.size());
326 
327  const double z = UgGridHelpers::cellCentroidCoordinate(G, *ci, nd-1);
328  p[c] = (z < split) ? f[up](z) : f[down](z);
329  }
330  }
331 
332  template <class Grid,
333  class Region,
334  class CellRange>
335  void
336  water(const Grid& G ,
337  const Region& reg ,
338  const std::array<double,2>& span ,
339  const double grav ,
340  double& po_woc,
341  const CellRange& cells ,
342  std::vector<double>& press )
343  {
344  using PhasePressODE::Water;
345  typedef Water<typename Region::CalcDensity> ODE;
346 
347  const PhaseUsage& pu = reg.phaseUsage();
348 
349  const int wix = PhaseIndex::water(pu);
350  const double T = 273.15 + 20; // standard temperature for now
351  ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
352 
353  double z0;
354  double p0;
355  if (reg.datum() > reg.zwoc()) {//Datum in water zone
356  z0 = reg.datum();
357  p0 = reg.pressure();
358  } else {
359  z0 = reg.zwoc();
360  p0 = po_woc - reg.pcow_woc(); // Water pressure at contact
361  }
362 
363  std::array<double,2> up = {{ z0, span[0] }};
364  std::array<double,2> down = {{ z0, span[1] }};
365 
366  typedef Details::RK4IVP<ODE> WPress;
367  std::array<WPress,2> wpress = {
368  {
369  WPress(drho, up , p0, 100)
370  ,
371  WPress(drho, down, p0, 100)
372  }
373  };
374 
375  assign(G, wpress, z0, cells, press);
376 
377  if (reg.datum() > reg.zwoc()) {
378  // Return oil pressure at contact
379  po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc();
380  }
381  }
382 
383  template <class Grid,
384  class Region,
385  class CellRange>
386  void
387  oil(const Grid& G ,
388  const Region& reg ,
389  const std::array<double,2>& span ,
390  const double grav ,
391  const CellRange& cells ,
392  std::vector<double>& press ,
393  double& po_woc,
394  double& po_goc)
395  {
396  using PhasePressODE::Oil;
397  typedef Oil<typename Region::CalcDensity,
398  typename Region::CalcDissolution> ODE;
399 
400  const PhaseUsage& pu = reg.phaseUsage();
401 
402  const int oix = PhaseIndex::oil(pu);
403  const int gix = PhaseIndex::gas(pu);
404  const double T = 273.15 + 20; // standard temperature for now
405  ODE drho(T,
406  reg.densityCalculator(),
407  reg.dissolutionCalculator(),
408  pu.num_phases, oix, gix, grav);
409 
410  double z0;
411  double p0;
412  if (reg.datum() > reg.zwoc()) {//Datum in water zone, po_woc given
413  z0 = reg.zwoc();
414  p0 = po_woc;
415  } else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, po_goc given
416  z0 = reg.zgoc();
417  p0 = po_goc;
418  } else { //Datum in oil zone
419  z0 = reg.datum();
420  p0 = reg.pressure();
421  }
422 
423  std::array<double,2> up = {{ z0, span[0] }};
424  std::array<double,2> down = {{ z0, span[1] }};
425 
426  typedef Details::RK4IVP<ODE> OPress;
427  std::array<OPress,2> opress = {
428  {
429  OPress(drho, up , p0, 100)
430  ,
431  OPress(drho, down, p0, 100)
432  }
433  };
434 
435  assign(G, opress, z0, cells, press);
436 
437  const double woc = reg.zwoc();
438  if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum
439  else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum
440  else { po_woc = p0; } // WOC *at* datum
441 
442  const double goc = reg.zgoc();
443  if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum
444  else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum
445  else { po_goc = p0; } // GOC *at* datum
446  }
447 
448  template <class Grid,
449  class Region,
450  class CellRange>
451  void
452  gas(const Grid& G ,
453  const Region& reg ,
454  const std::array<double,2>& span ,
455  const double grav ,
456  double& po_goc,
457  const CellRange& cells ,
458  std::vector<double>& press )
459  {
460  using PhasePressODE::Gas;
461  typedef Gas<typename Region::CalcDensity,
462  typename Region::CalcEvaporation> ODE;
463 
464  const PhaseUsage& pu = reg.phaseUsage();
465 
466  const int gix = PhaseIndex::gas(pu);
467  const int oix = PhaseIndex::oil(pu);
468 
469  const double T = 273.15 + 20; // standard temperature for now
470  ODE drho(T,
471  reg.densityCalculator(),
472  reg.evaporationCalculator(),
473  pu.num_phases, gix, oix, grav);
474 
475  double z0;
476  double p0;
477  if (reg.datum() < reg.zgoc()) {//Datum in gas zone
478  z0 = reg.datum();
479  p0 = reg.pressure();
480  } else {
481  z0 = reg.zgoc();
482  p0 = po_goc + reg.pcgo_goc(); // Gas pressure at contact
483  }
484 
485  std::array<double,2> up = {{ z0, span[0] }};
486  std::array<double,2> down = {{ z0, span[1] }};
487 
488  typedef Details::RK4IVP<ODE> GPress;
489  std::array<GPress,2> gpress = {
490  {
491  GPress(drho, up , p0, 100)
492  ,
493  GPress(drho, down, p0, 100)
494  }
495  };
496 
497  assign(G, gpress, z0, cells, press);
498 
499  if (reg.datum() < reg.zgoc()) {
500  // Return oil pressure at contact
501  po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc();
502  }
503  }
504  } // namespace PhasePressure
505 
506  template <class Grid,
507  class Region,
508  class CellRange>
509  void
510  equilibrateOWG(const Grid& G,
511  const Region& reg,
512  const double grav,
513  const std::array<double,2>& span,
514  const CellRange& cells,
515  std::vector< std::vector<double> >& press)
516  {
517  const PhaseUsage& pu = reg.phaseUsage();
518 
519  if (reg.datum() > reg.zwoc()) { // Datum in water zone
520  double po_woc = -1;
521  double po_goc = -1;
522 
523  if (PhaseUsed::water(pu)) {
524  const int wix = PhaseIndex::water(pu);
525  PhasePressure::water(G, reg, span, grav, po_woc,
526  cells, press[ wix ]);
527  }
528 
529  if (PhaseUsed::oil(pu)) {
530  const int oix = PhaseIndex::oil(pu);
531  PhasePressure::oil(G, reg, span, grav, cells,
532  press[ oix ], po_woc, po_goc);
533  }
534 
535  if (PhaseUsed::gas(pu)) {
536  const int gix = PhaseIndex::gas(pu);
537  PhasePressure::gas(G, reg, span, grav, po_goc,
538  cells, press[ gix ]);
539  }
540  } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone
541  double po_woc = -1;
542  double po_goc = -1;
543 
544  if (PhaseUsed::gas(pu)) {
545  const int gix = PhaseIndex::gas(pu);
546  PhasePressure::gas(G, reg, span, grav, po_goc,
547  cells, press[ gix ]);
548  }
549 
550  if (PhaseUsed::oil(pu)) {
551  const int oix = PhaseIndex::oil(pu);
552  PhasePressure::oil(G, reg, span, grav, cells,
553  press[ oix ], po_woc, po_goc);
554  }
555 
556  if (PhaseUsed::water(pu)) {
557  const int wix = PhaseIndex::water(pu);
558  PhasePressure::water(G, reg, span, grav, po_woc,
559  cells, press[ wix ]);
560  }
561  } else { // Datum in oil zone
562  double po_woc = -1;
563  double po_goc = -1;
564 
565  if (PhaseUsed::oil(pu)) {
566  const int oix = PhaseIndex::oil(pu);
567  PhasePressure::oil(G, reg, span, grav, cells,
568  press[ oix ], po_woc, po_goc);
569  }
570 
571  if (PhaseUsed::water(pu)) {
572  const int wix = PhaseIndex::water(pu);
573  PhasePressure::water(G, reg, span, grav, po_woc,
574  cells, press[ wix ]);
575  }
576 
577  if (PhaseUsed::gas(pu)) {
578  const int gix = PhaseIndex::gas(pu);
579  PhasePressure::gas(G, reg, span, grav, po_goc,
580  cells, press[ gix ]);
581  }
582  }
583  }
584  } // namespace Details
585 
586 
587  namespace Equil {
588 
589 
590  template <class Grid,
591  class Region,
592  class CellRange>
593  std::vector< std::vector<double> >
594  phasePressures(const Grid& G,
595  const Region& reg,
596  const CellRange& cells,
597  const double grav)
598  {
599  std::array<double,2> span =
600  {{ std::numeric_limits<double>::max() ,
601  -std::numeric_limits<double>::max() }}; // Symm. about 0.
602 
603  int ncell = 0;
604  {
605  // This code is only supported in three space dimensions
606  assert (UgGridHelpers::dimensions(G) == 3);
607 
608  const int nd = UgGridHelpers::dimensions(G);
609 
610  // Define vertical span as
611  //
612  // [minimum(node depth(cells)), maximum(node depth(cells))]
613  //
614  // Note: We use a sledgehammer approach--looping all
615  // the nodes of all the faces of all the 'cells'--to
616  // compute those bounds. This necessarily entails
617  // visiting some nodes (and faces) multiple times.
618  //
619  // Note: The implementation of 'RK4IVP<>' implicitly
620  // imposes the requirement that cell centroids are all
621  // within this vertical span. That requirement is not
622  // checked.
623  auto cell2Faces = UgGridHelpers::cell2Faces(G);
624  auto faceVertices = UgGridHelpers::face2Vertices(G);
625 
626  for (typename CellRange::const_iterator
627  ci = cells.begin(), ce = cells.end();
628  ci != ce; ++ci, ++ncell)
629  {
630  for (auto fi=cell2Faces[*ci].begin(),
631  fe=cell2Faces[*ci].end();
632  fi != fe;
633  ++fi)
634  {
635  for (auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end();
636  i != e; ++i)
637  {
638  const double z = UgGridHelpers::vertexCoordinates(G, *i)[nd-1];
639 
640  if (z < span[0]) { span[0] = z; }
641  if (z > span[1]) { span[1] = z; }
642  }
643  }
644  }
645  }
646  const int np = reg.phaseUsage().num_phases;
647 
648  typedef std::vector<double> pval;
649  std::vector<pval> press(np, pval(ncell, 0.0));
650 
651  const double zwoc = reg.zwoc ();
652  const double zgoc = reg.zgoc ();
653 
654  // make sure goc and woc is within the span for the phase pressure calculation
655  span[0] = std::min(span[0],zgoc);
656  span[1] = std::max(span[1],zwoc);
657 
658  Details::equilibrateOWG(G, reg, grav, span, cells, press);
659 
660  return press;
661  }
662 
663  template <class Grid,
664  class Region,
665  class CellRange>
666  std::vector<double>
667  temperature(const Grid& /* G */,
668  const Region& /* reg */,
669  const CellRange& cells)
670  {
671  // use the standard temperature for everything for now
672  return std::vector<double>(cells.size(), 273.15 + 20.0);
673  }
674 
675  template <class Grid, class Region, class CellRange>
676  std::vector< std::vector<double> >
677  phaseSaturations(const Grid& G,
678  const Region& reg,
679  const CellRange& cells,
680  BlackoilPropertiesInterface& props,
681  const std::vector<double> swat_init,
682  std::vector< std::vector<double> >& phase_pressures)
683  {
684  if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
685  OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
686  }
687 
688  std::vector< std::vector<double> > phase_saturations = phase_pressures; // Just to get the right size.
689  double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
690  double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
691 
692  const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
693  const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
694  const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
695  const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
696  const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
697  std::vector<double>::size_type local_index = 0;
698  for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
699  const int cell = *ci;
700  props.satRange(1, &cell, smin, smax);
701  // Find saturations from pressure differences by
702  // inverting capillary pressure functions.
703  double sw = 0.0;
704  if (water) {
705  if (isConstPc(props,waterpos,cell)){
706  const int nd = UgGridHelpers::dimensions(G);
707  const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G,
708  cell,
709  nd-1);
710  sw = satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,false);
711  phase_saturations[waterpos][local_index] = sw;
712  }
713  else{
714  const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
715  if (swat_init.empty()) { // Invert Pc to find sw
716  sw = satFromPc(props, waterpos, cell, pcov);
717  phase_saturations[waterpos][local_index] = sw;
718  } else { // Scale Pc to reflect imposed sw
719  sw = swat_init[cell];
720  props.swatInitScaling(cell, pcov, sw);
721  phase_saturations[waterpos][local_index] = sw;
722  }
723  }
724  }
725  double sg = 0.0;
726  if (gas) {
727  if (isConstPc(props,gaspos,cell)){
728  const int nd = UgGridHelpers::dimensions(G);
729  const double cellDepth = UgGridHelpers::cellCentroidCoordinate(G,
730  cell,
731  nd-1);
732  sg = satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,true);
733  phase_saturations[gaspos][local_index] = sg;
734  }
735  else{
736  // Note that pcog is defined to be (pg - po), not (po - pg).
737  const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
738  const double increasing = true; // pcog(sg) expected to be increasing function
739  sg = satFromPc(props, gaspos, cell, pcog, increasing);
740  phase_saturations[gaspos][local_index] = sg;
741  }
742  }
743  if (gas && water && (sg + sw > 1.0)) {
744  // Overlapping gas-oil and oil-water transition
745  // zones can lead to unphysical saturations when
746  // treated as above. Must recalculate using gas-water
747  // capillary pressure.
748  const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
749  if (! swat_init.empty()) {
750  // Re-scale Pc to reflect imposed sw for vanishing oil phase.
751  // This seems consistent with ecl, and fails to honour
752  // swat_init in case of non-trivial gas-oil cap pressure.
753  props.swatInitScaling(cell, pcgw, sw);
754  }
755  sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
756  sg = 1.0 - sw;
757  phase_saturations[waterpos][local_index] = sw;
758  phase_saturations[gaspos][local_index] = sg;
759  // Adjust oil pressure according to gas saturation and cap pressure
760  double pc[BlackoilPhases::MaxNumPhases];
761  double sat[BlackoilPhases::MaxNumPhases];
762  sat[waterpos] = sw;
763  sat[gaspos] = sg;
764  sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
765  props.capPress(1, sat, &cell, pc, 0);
766  phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
767  }
768  phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
769 
770  // Adjust phase pressures for max and min saturation ...
771  double pc[BlackoilPhases::MaxNumPhases];
772  double sat[BlackoilPhases::MaxNumPhases];
773  double threshold_sat = 1.0e-6;
774 
775  sat[waterpos] = smax[waterpos];
776  sat[gaspos] = smax[gaspos];
777  sat[oilpos] = 1.0 - sat[waterpos] - sat[gaspos];
778  if (sw > smax[waterpos]-threshold_sat ) {
779  sat[waterpos] = smax[waterpos];
780  props.capPress(1, sat, &cell, pc, 0);
781  phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
782  } else if (sg > smax[gaspos]-threshold_sat) {
783  sat[gaspos] = smax[gaspos];
784  props.capPress(1, sat, &cell, pc, 0);
785  phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
786  }
787  if (sg < smin[gaspos]+threshold_sat) {
788  sat[gaspos] = smin[gaspos];
789  props.capPress(1, sat, &cell, pc, 0);
790  phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
791  }
792  if (sw < smin[waterpos]+threshold_sat) {
793  sat[waterpos] = smin[waterpos];
794  props.capPress(1, sat, &cell, pc, 0);
795  phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
796  }
797  }
798  return phase_saturations;
799  }
800 
801 
820  template <class Grid, class CellRangeType>
821  std::vector<double> computeRs(const Grid& grid,
822  const CellRangeType& cells,
823  const std::vector<double> oil_pressure,
824  const std::vector<double>& temperature,
825  const Miscibility::RsFunction& rs_func,
826  const std::vector<double> gas_saturation)
827  {
828  assert(UgGridHelpers::dimensions(grid) == 3);
829  std::vector<double> rs(cells.size());
830  int count = 0;
831  for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
832  const double depth = UgGridHelpers::cellCentroidCoordinate(grid, *it, 2);
833  rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
834  }
835  return rs;
836  }
837 
838  } // namespace Equil
839 
840 
841  namespace Details
842  {
846  inline std::vector<double>
847  convertSats(const std::vector< std::vector<double> >& sat)
848  {
849  const auto np = sat.size();
850  const auto nc = sat[0].size();
851 
852  std::vector<double> s(np * nc);
853 
854  for (decltype(sat.size()) p = 0; p < np; ++p) {
855  const auto& sat_p = sat[p];
856  double* sp = & s[0*nc + p];
857 
858  for (decltype(sat[0].size()) c = 0;
859  c < nc; ++c, sp += np)
860  {
861  *sp = sat_p[c];
862  }
863  }
864 
865  return s;
866  }
867  } // namespace Details
868 
869 
885  template<class Grid>
886  void initStateEquil(const Grid& grid,
887  BlackoilPropertiesInterface& props,
888  const Opm::DeckConstPtr deck,
889  const Opm::EclipseStateConstPtr eclipseState,
890  const double gravity,
891  BlackoilState& state)
892  {
893  typedef Equil::DeckDependent::InitialStateComputer ISC;
894  ISC isc(props, deck, eclipseState, grid, gravity);
895  const auto pu = props.phaseUsage();
896  const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
897  ? pu.phase_pos[BlackoilPhases::Liquid]
898  : pu.phase_pos[BlackoilPhases::Aqua];
899  state.pressure() = isc.press()[ref_phase];
900  state.saturation() = Details::convertSats(isc.saturation());
901  state.gasoilratio() = isc.rs();
902  state.rv() = isc.rv();
903  initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
904  }
905 
906 
907 
908 } // namespace Opm
909 
910 #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
Opm::Details::PhasePressODE::Water
Definition: initStateEquil_impl.hpp:110
Opm::Details::PhaseIndex::oil
int oil(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:283
Opm::Details::PhaseIndex::gas
int gas(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:294
Opm::parameter::split
std::pair< std::string, std::string > split(const std::string &name)
Opm::Details::PhasePressODE::Water::Water
Water(const double temp, const Density &rho, const int np, const int ix, const double norm_grav)
Definition: initStateEquil_impl.hpp:112
Opm::initStateEquil
void initStateEquil(const Grid &grid, const BlackoilPropertiesInterface &props, const Opm::DeckConstPtr deck, const Opm::EclipseStateConstPtr eclipseState, const double gravity, BlackoilState &state)
Opm::Details::PhasePressure::gas
void gas(const Grid &G, const Region &reg, const std::array< double, 2 > &span, const double grav, double &po_goc, const CellRange &cells, std::vector< double > &press)
Definition: initStateEquil_impl.hpp:452
Opm::Details::PhasePressODE::Oil::operator()
double operator()(const double depth, const double press) const
Definition: initStateEquil_impl.hpp:171
Opm::Details::equilibrateOWG
void equilibrateOWG(const Grid &G, const Region &reg, const double grav, const std::array< double, 2 > &span, const CellRange &cells, std::vector< std::vector< double > > &press)
Definition: initStateEquil_impl.hpp:510
Opm
Definition: AnisotropicEikonal.hpp:43
Opm::Details::PhasePressODE::Oil::Oil
Oil(const double temp, const Density &rho, const RS &rs, const int np, const int oix, const int gix, const double norm_grav)
Definition: initStateEquil_impl.hpp:152
Opm::BlackoilPhases::Vapour
Definition: BlackoilPhases.hpp:32
Opm::BlackoilState::rv
std::vector< double > & rv()
Definition: BlackoilState.hpp:55
Opm::Details::PhaseUsed::gas
bool gas(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:264
Opm::SimulatorState::saturation
std::vector< double > & saturation()
Definition: SimulatorState.hpp:60
Opm::Details::PhasePressure::assign
void assign(const Grid &G, const std::array< PressFunction, 2 > &f, const double split, const CellRange &cells, std::vector< double > &p)
Definition: initStateEquil_impl.hpp:310
Opm::Details::PhasePressODE::Oil
Definition: initStateEquil_impl.hpp:150
Opm::UgGridHelpers::face2Vertices
Face2VerticesTraits< UnstructuredGrid >::Type face2Vertices(const UnstructuredGrid &grid)
Get the face to vertices mapping of a grid.
Opm::BlackoilPropertiesInterface::capPress
virtual void capPress(const int n, const double *s, const int *cells, double *pc, double *dpcds) const =0
Opm::PhaseUsage::phase_used
int phase_used[MaxNumPhases]
Definition: BlackoilPhases.hpp:39
Opm::Equil::satFromDepth
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
Opm::Details::PhasePressure::water
void water(const Grid &G, const Region &reg, const std::array< double, 2 > &span, const double grav, double &po_woc, const CellRange &cells, std::vector< double > &press)
Definition: initStateEquil_impl.hpp:336
Opm::Equil::isConstPc
bool isConstPc(const BlackoilPropertiesInterface &props, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers.hpp:890
Opm::Details::PhaseUsed::water
bool water(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:252
Opm::UgGridHelpers::numCells
int numCells(const UnstructuredGrid &grid)
Get the number of cells of a grid.
Opm::Details::PhasePressODE::Gas
Definition: initStateEquil_impl.hpp:200
Opm::BlackoilPropertiesInterface
Definition: BlackoilPropertiesInterface.hpp:37
Opm::Details::PhasePressure::oil
void oil(const Grid &G, const Region &reg, const std::array< double, 2 > &span, const double grav, const CellRange &cells, std::vector< double > &press, double &po_woc, double &po_goc)
Definition: initStateEquil_impl.hpp:387
Opm::PhaseUsage::phase_pos
int phase_pos[MaxNumPhases]
Definition: BlackoilPhases.hpp:40
Opm::Details::PhasePressODE::Gas::Gas
Gas(const double temp, const Density &rho, const RV &rv, const int np, const int gix, const int oix, const double norm_grav)
Definition: initStateEquil_impl.hpp:202
Opm::Equil::phasePressures
std::vector< std::vector< double > > phasePressures(const Grid &G, const Region &reg, const CellRange &cells, const double grav)
Definition: initStateEquil_impl.hpp:594
Opm::BlackoilPropertiesInterface::satRange
virtual void satRange(const int n, const int *cells, double *smin, double *smax) const =0
Opm::PhaseUsage::num_phases
int num_phases
Definition: BlackoilPhases.hpp:38
Opm::BlackoilState
Simulator state for a blackoil simulator.
Definition: BlackoilState.hpp:33
Opm::UgGridHelpers::cellCentroidCoordinate
double cellCentroidCoordinate(const UnstructuredGrid &grid, int cell_index, int coordinate)
Get a coordinate of a specific cell centroid.
Opm::UgGridHelpers::dimensions
int dimensions(const UnstructuredGrid &grid)
Get the dimensions of a grid.
Opm::Details::convertSats
std::vector< double > convertSats(const std::vector< std::vector< double > > &sat)
Definition: initStateEquil_impl.hpp:847
Opm::BlackoilPropertiesInterface::phaseUsage
virtual PhaseUsage phaseUsage() const =0
Opm::Equil::satFromPc
double satFromPc(const BlackoilPropertiesInterface &props, const int phase, const int cell, const double target_pc, const bool increasing=false)
Definition: EquilibrationHelpers.hpp:760
Opm::unit::gravity
const double gravity
Definition: Units.hpp:120
Opm::Details::RK4IVP::operator()
double operator()(const double x) const
Definition: initStateEquil_impl.hpp:75
BlackoilPhases.hpp
Opm::Equil::computeRs
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)
Definition: initStateEquil_impl.hpp:821
GridHelpers.hpp
Opm::BlackoilPhases::Liquid
Definition: BlackoilPhases.hpp:32
Opm::Details::RK4IVP
Definition: initStateEquil_impl.hpp:39
Opm::Details::PhasePressODE::Gas::operator()
double operator()(const double depth, const double press) const
Definition: initStateEquil_impl.hpp:221
grid
const UnstructuredGrid & grid
Definition: ColumnExtract.hpp:31
Opm::BlackoilPhases::MaxNumPhases
static const int MaxNumPhases
Definition: BlackoilPhases.hpp:30
Opm::Details::PhaseUsed::oil
bool oil(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:258
Opm::UgGridHelpers::cell2Faces
Cell2FacesTraits< UnstructuredGrid >::Type cell2Faces(const UnstructuredGrid &grid)
Get the cell to faces mapping of a grid.
Opm::BlackoilState::gasoilratio
std::vector< double > & gasoilratio()
Definition: BlackoilState.hpp:54
Opm::initBlackoilSurfvolUsingRSorRV
void initBlackoilSurfvolUsingRSorRV(const UnstructuredGrid &grid, const Props &props, State &state)
Definition: initState_impl.hpp:724
Opm::Details::PhaseIndex::water
int water(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:272
Opm::BlackoilPropertiesInterface::swatInitScaling
virtual void swatInitScaling(const int cell, const double pcow, double &swat)=0
Opm::Details::RK4IVP::RK4IVP
RK4IVP(const RHS &f, const std::array< double, 2 > &span, const double y0, const int N)
Definition: initStateEquil_impl.hpp:41
Opm::Details::PhasePressODE::Water::operator()
double operator()(const double, const double press) const
Definition: initStateEquil_impl.hpp:127
Opm::SimulatorState::pressure
std::vector< double > & pressure()
Definition: SimulatorState.hpp:56
Opm::UgGridHelpers::vertexCoordinates
const double * vertexCoordinates(const UnstructuredGrid &grid, int index)
Get the coordinates of a vertex of the grid.
Opm::BlackoilPhases::Aqua
Definition: BlackoilPhases.hpp:32
Opm::Equil::temperature
std::vector< double > temperature(const Grid &, const Region &, const CellRange &cells)
Definition: initStateEquil_impl.hpp:667
Opm::PhaseUsage
Definition: BlackoilPhases.hpp:36
Opm::Equil::Miscibility::RsFunction
Definition: EquilibrationHelpers.hpp:165
Opm::Equil::satFromSumOfPcs
double satFromSumOfPcs(const BlackoilPropertiesInterface &props, const int phase1, const int phase2, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:835
Opm::Equil::phaseSaturations
std::vector< std::vector< double > > phaseSaturations(const Grid &G, const Region &reg, const CellRange &cells, BlackoilPropertiesInterface &props, const std::vector< double > swat_init, std::vector< std::vector< double > > &phase_pressures)
Definition: initStateEquil_impl.hpp:677