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