23#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
24#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
26#include <opm/core/grid.h>
27#include <opm/core/grid/GridHelpers.hpp>
29#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
43 const std::array<double,2>& span,
49 const double h = stepsize();
50 const double h2 = h / 2;
51 const double h6 = h / 6;
57 f_.push_back(f(span_[0], y0));
59 for (
int i = 0; i < N; ++i) {
60 const double x = span_[0] + i*h;
61 const double y = y_.back();
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);
68 y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
69 f_.push_back(f(x + h, y_.back()));
72 assert (y_.size() == std::vector<double>::size_type(N + 1));
80 const double h = stepsize();
81 int i = (x - span_[0]) / h;
82 const double t = (x - (span_[0] + i*h)) / h;
86 if (N_ <= i) { i = N_ - 1; }
88 const double y0 = y_[i], y1 = y_[i + 1];
89 const double f0 = f_[i], f1 = f_[i + 1];
91 double u = (1 - 2*t) * (y1 - y0);
92 u += h * ((t - 1)*f0 + t*f1);
94 u += (1 - t)*y0 + t*y1;
101 std::array<double,2> span_;
102 std::vector<double> y_;
103 std::vector<double> f_;
106 stepsize()
const {
return (span_[1] - span_[0]) / N_; }
109 namespace PhasePressODE {
110 template <
class Flu
idSystem>
114 const int pvtRegionIdx,
115 const double norm_grav)
117 , pvtRegionIdx_(pvtRegionIdx)
124 const double press)
const
126 return this->density(press) * g_;
131 const int pvtRegionIdx_;
135 density(
const double press)
const
137 double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
138 rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_);
143 template <
class Flu
idSystem,
class RS>
148 const int pvtRegionIdx,
149 const double norm_grav)
152 , pvtRegionIdx_(pvtRegionIdx)
159 const double press)
const
161 return this->density(depth, press) * g_;
167 const int pvtRegionIdx_;
171 density(
const double depth,
172 const double press)
const
174 double rs = rs_(depth, press, temp_);
176 if ( !FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press) ) {
177 bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
179 bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs);
181 double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
182 if (FluidSystem::enableDissolvedGas()) {
183 rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
190 template <
class Flu
idSystem,
class RV>
195 const int pvtRegionIdx,
196 const double norm_grav)
199 , pvtRegionIdx_(pvtRegionIdx)
206 const double press)
const
208 return this->density(depth, press) * g_;
214 const int pvtRegionIdx_;
218 density(
const double depth,
219 const double press)
const
221 double rv = rv_(depth, press, temp_);
223 if ( !FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) ) {
224 bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press);
226 bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv);
228 double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_);
229 if (FluidSystem::enableVaporizedOil()) {
230 rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_);
239 namespace PhasePressure {
240 template <
class Grid,
245 const std::array<PressFunction, 2>& f ,
247 const CellRange& cells,
248 std::vector<double>& p )
251 enum { up = 0, down = 1 };
253 std::vector<double>::size_type c = 0;
254 for (
typename CellRange::const_iterator
255 ci = cells.begin(), ce = cells.end();
258 assert (c < p.size());
260 const double z = UgGridHelpers::cellCenterDepth(G, *ci);
261 p[c] = (z < split) ? f[up](z) : f[down](z);
265 template <
class FluidSystem,
272 const std::array<double,2>& span ,
275 const CellRange& cells ,
276 std::vector<double>& press )
279 typedef Water<FluidSystem> ODE;
281 const double T = 273.15 + 20;
282 ODE drho(T, reg.pvtIdx() , grav);
286 if (reg.datum() > reg.zwoc()) {
291 p0 = po_woc - reg.pcow_woc();
294 std::array<double,2> up = {{ z0, span[0] }};
295 std::array<double,2> down = {{ z0, span[1] }};
298 std::array<WPress,2> wpress = {
300 WPress(drho, up , p0, 2000)
302 WPress(drho, down, p0, 2000)
306 assign(G, wpress, z0, cells, press);
308 if (reg.datum() > reg.zwoc()) {
310 po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc();
314 template <
class FluidSystem,
321 const std::array<double,2>& span ,
323 const CellRange& cells ,
324 std::vector<double>& press ,
329 typedef Oil<FluidSystem, typename Region::CalcDissolution> ODE;
331 const double T = 273.15 + 20;
332 ODE drho(T, reg.dissolutionCalculator(),
337 if (reg.datum() > reg.zwoc()) {
340 }
else if (reg.datum() < reg.zgoc()) {
348 std::array<double,2> up = {{ z0, span[0] }};
349 std::array<double,2> down = {{ z0, span[1] }};
352 std::array<OPress,2> opress = {
354 OPress(drho, up , p0, 2000)
356 OPress(drho, down, p0, 2000)
360 assign(G, opress, z0, cells, press);
362 const double woc = reg.zwoc();
363 if (z0 > woc) { po_woc = opress[0](woc); }
364 else if (z0 < woc) { po_woc = opress[1](woc); }
365 else { po_woc = p0; }
367 const double goc = reg.zgoc();
368 if (z0 > goc) { po_goc = opress[0](goc); }
369 else if (z0 < goc) { po_goc = opress[1](goc); }
370 else { po_goc = p0; }
373 template <
class FluidSystem,
380 const std::array<double,2>& span ,
383 const CellRange& cells ,
384 std::vector<double>& press )
387 typedef Gas<FluidSystem, typename Region::CalcEvaporation> ODE;
389 const double T = 273.15 + 20;
390 ODE drho(T, reg.evaporationCalculator(),
395 if (reg.datum() < reg.zgoc()) {
400 p0 = po_goc + reg.pcgo_goc();
403 std::array<double,2> up = {{ z0, span[0] }};
404 std::array<double,2> down = {{ z0, span[1] }};
407 std::array<GPress,2> gpress = {
409 GPress(drho, up , p0, 2000)
411 GPress(drho, down, p0, 2000)
415 assign(G, gpress, z0, cells, press);
417 if (reg.datum() < reg.zgoc()) {
419 po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc();
424 template <
class FluidSystem,
432 const std::array<double,2>& span,
433 const CellRange& cells,
434 std::vector< std::vector<double> >& press)
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;
443 if (reg.datum() > reg.zwoc()) {
448 PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
449 cells, press[ waterpos ]);
453 PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
454 press[ oilpos ], po_woc, po_goc);
458 PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
459 cells, press[ gaspos ]);
461 }
else if (reg.datum() < reg.zgoc()) {
466 PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
467 cells, press[ gaspos ]);
471 PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
472 press[ oilpos ], po_woc, po_goc);
476 PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
477 cells, press[ waterpos ]);
484 PhasePressure::oil<FluidSystem>(G, reg, span, grav, cells,
485 press[ oilpos ], po_woc, po_goc);
489 PhasePressure::water<FluidSystem>(G, reg, span, grav, po_woc,
490 cells, press[ waterpos ]);
494 PhasePressure::gas<FluidSystem>(G, reg, span, grav, po_goc,
495 cells, press[ gaspos ]);
505 template <
class FluidSystem,
509 std::vector< std::vector<double> >
512 const CellRange& cells,
515 std::array<double,2> span =
516 {{ std::numeric_limits<double>::max() ,
517 -std::numeric_limits<double>::max() }};
522 assert (UgGridHelpers::dimensions(G) == 3);
524 const int nd = UgGridHelpers::dimensions(G);
539 auto cell2Faces = UgGridHelpers::cell2Faces(G);
540 auto faceVertices = UgGridHelpers::face2Vertices(G);
542 for (
typename CellRange::const_iterator
543 ci = cells.begin(), ce = cells.end();
544 ci != ce; ++ci, ++ncell)
546 for (
auto fi=cell2Faces[*ci].begin(),
547 fe=cell2Faces[*ci].end();
551 for (
auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end();
554 const double z = UgGridHelpers::vertexCoordinates(G, *i)[nd-1];
556 if (z < span[0]) { span[0] = z; }
557 if (z > span[1]) { span[1] = z; }
562 const int np = FluidSystem::numPhases;
564 typedef std::vector<double> pval;
565 std::vector<pval> press(np, pval(ncell, 0.0));
567 const double zwoc = reg.zwoc ();
568 const double zgoc = reg.zgoc ();
571 span[0] = std::min(span[0],zgoc);
572 span[1] = std::max(span[1],zwoc);
574 Details::equilibrateOWG<FluidSystem>(G, reg, grav, span, cells, press);
579 template <
class Grid,
585 const CellRange& cells)
588 return std::vector<double>(cells.size(), 273.15 + 20.0);
591 template <
class Flu
idSystem,
class Gr
id,
class Region,
class CellRange,
class MaterialLawManager>
592 std::vector< std::vector<double> >
595 const CellRange& cells,
596 MaterialLawManager& materialLawManager,
597 const std::vector<double> swat_init,
598 std::vector< std::vector<double> >& phase_pressures)
600 if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
601 OPM_THROW(std::runtime_error,
"Cannot initialise: not handling water-gas cases.");
604 std::vector< std::vector<double> > phase_saturations = phase_pressures;
607 typedef Opm::SimpleModularFluidState<double,
621 typedef typename MaterialLawManager::MaterialLaw MaterialLaw;
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);
639 if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::waterPhaseIdx, cell)){
640 const double cellDepth = UgGridHelpers::cellCenterDepth(G,
642 sw = satFromDepth<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,cellDepth,reg.zwoc(),waterpos,cell,
false);
643 phase_saturations[waterpos][local_index] = sw;
646 const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
647 if (swat_init.empty()) {
648 sw = satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, waterpos, cell, pcov);
649 phase_saturations[waterpos][local_index] = sw;
651 sw = swat_init[cell];
652 sw = materialLawManager.applySwatinit(cell, pcov, sw);
653 phase_saturations[waterpos][local_index] = sw;
659 if (isConstPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,FluidSystem::gasPhaseIdx,cell)){
660 const double cellDepth = UgGridHelpers::cellCenterDepth(G,
662 sg = satFromDepth<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager,cellDepth,reg.zgoc(),gaspos,cell,
true);
663 phase_saturations[gaspos][local_index] = sg;
667 const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
668 const double increasing =
true;
669 sg = satFromPc<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, gaspos, cell, pcog, increasing);
670 phase_saturations[gaspos][local_index] = sg;
673 if (
gas &&
water && (sg + sw > 1.0)) {
678 const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
679 if (! swat_init.empty()) {
683 sw = materialLawManager.applySwatinit(cell, pcgw, sw);
685 sw = satFromSumOfPcs<FluidSystem, MaterialLaw, MaterialLawManager>(materialLawManager, waterpos, gaspos, cell, pcgw);
687 phase_saturations[waterpos][local_index] = sw;
688 phase_saturations[gaspos][local_index] = sg;
690 fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw);
693 fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
695 fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - sw - sg);
696 fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg);
698 double pC[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;
703 phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
706 double threshold_sat = 1.0e-6;
709 double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 };
711 double swu = scaledDrainageInfo.Swu;
712 fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu);
716 double sgu = scaledDrainageInfo.Sgu;
717 fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu);
720 fluidState.setSaturation(FluidSystem::oilPhaseIdx, so);
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;
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;
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;
746 return phase_saturations;
768 template <
class Gr
id,
class CellRangeType>
770 const CellRangeType& cells,
771 const std::vector<double> oil_pressure,
774 const std::vector<double> gas_saturation)
776 assert(UgGridHelpers::dimensions(grid) == 3);
777 std::vector<double> rs(cells.size());
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]);
Definition: initStateEquil_impl.hpp:191
double operator()(const double depth, const double press) const
Definition: initStateEquil_impl.hpp:205
Gas(const double temp, const RV &rv, const int pvtRegionIdx, const double norm_grav)
Definition: initStateEquil_impl.hpp:193
Definition: initStateEquil_impl.hpp:144
double operator()(const double depth, const double press) const
Definition: initStateEquil_impl.hpp:158
Oil(const double temp, const RS &rs, const int pvtRegionIdx, const double norm_grav)
Definition: initStateEquil_impl.hpp:146
Definition: initStateEquil_impl.hpp:111
Water(const double temp, const int pvtRegionIdx, const double norm_grav)
Definition: initStateEquil_impl.hpp:113
double operator()(const double, const double press) const
Definition: initStateEquil_impl.hpp:123
Definition: initStateEquil_impl.hpp:40
RK4IVP(const RHS &f, const std::array< double, 2 > &span, const double y0, const int N)
Definition: initStateEquil_impl.hpp:42
double operator()(const double x) const
Definition: initStateEquil_impl.hpp:76
Definition: EquilibrationHelpers.hpp:116
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
void gas(const Grid &G, const Region ®, const std::array< double, 2 > &span, const double grav, double &po_goc, const CellRange &cells, std::vector< double > &press)
Definition: initStateEquil_impl.hpp:378
void water(const Grid &G, const Region ®, const std::array< double, 2 > &span, const double grav, double &po_woc, const CellRange &cells, std::vector< double > &press)
Definition: initStateEquil_impl.hpp:270
void oil(const Grid &G, const Region ®, 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
void equilibrateOWG(const Grid &G, const Region ®, const double grav, const std::array< double, 2 > &span, const CellRange &cells, std::vector< std::vector< double > > &press)
Definition: initStateEquil_impl.hpp:429
Opm::SimpleModularFluidState< double, 3, 3, FluidSystemSimple, false, false, false, false, true, false, false, false > SatOnlyFluidState
Definition: EquilibrationHelpers.hpp:104
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
std::vector< std::vector< double > > phaseSaturations(const Grid &G, const Region ®, const CellRange &cells, MaterialLawManager &materialLawManager, const std::vector< double > swat_init, std::vector< std::vector< double > > &phase_pressures)
Definition: initStateEquil_impl.hpp:593
std::vector< std::vector< double > > phasePressures(const Grid &G, const Region ®, const CellRange &cells, const double grav)
Definition: initStateEquil_impl.hpp:510
std::vector< double > temperature(const Grid &, const Region &, const CellRange &cells)
Definition: initStateEquil_impl.hpp:583
Definition: AnisotropicEikonal.hpp:44