22 #ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
23 #define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
25 #include <opm/core/grid.h>
26 #include <opm/core/grid/GridHelpers.hpp>
27 #include <opm/core/props/BlackoilPhases.hpp>
42 const std::array<double,2>& span,
48 const double h = stepsize();
49 const double h2 = h / 2;
50 const double h6 = h / 6;
56 f_.push_back(f(span_[0], y0));
58 for (
int i = 0; i < N; ++i) {
59 const double x = span_[0] + i*h;
60 const double y = y_.back();
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);
67 y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
68 f_.push_back(f(x + h, y_.back()));
71 assert (y_.size() == std::vector<double>::size_type(N + 1));
75 operator()(
const double x)
const
79 const double h = stepsize();
80 int i = (x - span_[0]) / h;
81 const double t = (x - (span_[0] + i*h)) / h;
85 if (N_ <= i) { i = N_ - 1; }
87 const double y0 = y_[i], y1 = y_[i + 1];
88 const double f0 = f_[i], f1 = f_[i + 1];
90 double u = (1 - 2*t) * (y1 - y0);
91 u += h * ((t - 1)*f0 + t*f1);
93 u += (1 - t)*y0 + t*y1;
100 std::array<double,2> span_;
101 std::vector<double> y_;
102 std::vector<double> f_;
105 stepsize()
const {
return (span_[1] - span_[0]) / N_; }
108 namespace PhasePressODE {
109 template <
class Density>
112 Water(
const double temp,
116 const double norm_grav)
127 operator()(
const double ,
128 const double press)
const
130 return this->density(press) * g_;
136 std::vector<double> svol_;
141 density(
const double press)
const
143 const std::vector<double>& rho = rho_(press, temp_, svol_);
149 template <
class Density,
class RS>
152 Oil(
const double temp,
158 const double norm_grav)
171 operator()(
const double depth,
172 const double press)
const
174 return this->density(depth, press) * g_;
181 mutable std::vector<double> svol_;
187 density(
const double depth,
188 const double press)
const
191 svol_[gix_] = rs_(depth, press, temp_);
194 const std::vector<double>& rho = rho_(press, temp_, svol_);
199 template <
class Density,
class RV>
202 Gas(
const double temp,
208 const double norm_grav)
221 operator()(
const double depth,
222 const double press)
const
224 return this->density(depth, press) * g_;
231 mutable std::vector<double> svol_;
237 density(
const double depth,
238 const double press)
const
241 svol_[oix_] = rv_(depth, press, temp_);
244 const std::vector<double>& rho = rho_(press, temp_, svol_);
250 namespace PhaseUsed {
254 return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
260 return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
264 gas(
const PhaseUsage& pu)
266 return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
270 namespace PhaseIndex {
272 water(
const PhaseUsage& pu)
275 if (PhaseUsed::water(pu)) {
276 i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
283 oil(
const PhaseUsage& pu)
286 if (PhaseUsed::oil(pu)) {
287 i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
294 gas(
const PhaseUsage& pu)
297 if (PhaseUsed::gas(pu)) {
298 i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
305 namespace PhasePressure {
306 template <
class Grid,
310 assign(
const Grid& G ,
311 const std::array<PressFunction, 2>& f ,
313 const CellRange& cells,
314 std::vector<double>& p )
317 enum { up = 0, down = 1 };
319 std::vector<double>::size_type c = 0;
320 for (
typename CellRange::const_iterator
321 ci = cells.begin(), ce = cells.end();
324 assert (c < p.size());
326 const double z = UgGridHelpers::cellCenterDepth(G, *ci);
327 p[c] = (z < split) ? f[up](z) : f[down](z);
331 template <
class Grid,
335 water(
const Grid& G ,
337 const std::array<double,2>& span ,
340 const CellRange& cells ,
341 std::vector<double>& press )
343 using PhasePressODE::Water;
344 typedef Water<typename Region::CalcDensity> ODE;
346 const PhaseUsage& pu = reg.phaseUsage();
348 const int wix = PhaseIndex::water(pu);
349 const double T = 273.15 + 20;
350 ODE drho(T, reg.densityCalculator(), pu.num_phases, wix, grav);
354 if (reg.datum() > reg.zwoc()) {
359 p0 = po_woc - reg.pcow_woc();
362 std::array<double,2> up = {{ z0, span[0] }};
363 std::array<double,2> down = {{ z0, span[1] }};
365 typedef Details::RK4IVP<ODE> WPress;
366 std::array<WPress,2> wpress = {
368 WPress(drho, up , p0, 2000)
370 WPress(drho, down, p0, 2000)
374 assign(G, wpress, z0, cells, press);
376 if (reg.datum() > reg.zwoc()) {
378 po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc();
382 template <
class Grid,
388 const std::array<double,2>& span ,
390 const CellRange& cells ,
391 std::vector<double>& press ,
395 using PhasePressODE::Oil;
396 typedef Oil<
typename Region::CalcDensity,
397 typename Region::CalcDissolution> ODE;
399 const PhaseUsage& pu = reg.phaseUsage();
401 const int oix = PhaseIndex::oil(pu);
402 const int gix = PhaseIndex::gas(pu);
403 const double T = 273.15 + 20;
405 reg.densityCalculator(),
406 reg.dissolutionCalculator(),
407 pu.num_phases, oix, gix, grav);
411 if (reg.datum() > reg.zwoc()) {
414 }
else if (reg.datum() < reg.zgoc()) {
422 std::array<double,2> up = {{ z0, span[0] }};
423 std::array<double,2> down = {{ z0, span[1] }};
425 typedef Details::RK4IVP<ODE> OPress;
426 std::array<OPress,2> opress = {
428 OPress(drho, up , p0, 2000)
430 OPress(drho, down, p0, 2000)
434 assign(G, opress, z0, cells, press);
436 const double woc = reg.zwoc();
437 if (z0 > woc) { po_woc = opress[0](woc); }
438 else if (z0 < woc) { po_woc = opress[1](woc); }
439 else { po_woc = p0; }
441 const double goc = reg.zgoc();
442 if (z0 > goc) { po_goc = opress[0](goc); }
443 else if (z0 < goc) { po_goc = opress[1](goc); }
444 else { po_goc = p0; }
447 template <
class Grid,
453 const std::array<double,2>& span ,
456 const CellRange& cells ,
457 std::vector<double>& press )
459 using PhasePressODE::Gas;
460 typedef Gas<
typename Region::CalcDensity,
461 typename Region::CalcEvaporation> ODE;
463 const PhaseUsage& pu = reg.phaseUsage();
465 const int gix = PhaseIndex::gas(pu);
466 const int oix = PhaseIndex::oil(pu);
468 const double T = 273.15 + 20;
470 reg.densityCalculator(),
471 reg.evaporationCalculator(),
472 pu.num_phases, gix, oix, grav);
476 if (reg.datum() < reg.zgoc()) {
481 p0 = po_goc + reg.pcgo_goc();
484 std::array<double,2> up = {{ z0, span[0] }};
485 std::array<double,2> down = {{ z0, span[1] }};
487 typedef Details::RK4IVP<ODE> GPress;
488 std::array<GPress,2> gpress = {
490 GPress(drho, up , p0, 2000)
492 GPress(drho, down, p0, 2000)
496 assign(G, gpress, z0, cells, press);
498 if (reg.datum() < reg.zgoc()) {
500 po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc();
505 template <
class Grid,
509 equilibrateOWG(
const Grid& G,
512 const std::array<double,2>& span,
513 const CellRange& cells,
514 std::vector< std::vector<double> >& press)
516 const PhaseUsage& pu = reg.phaseUsage();
518 if (reg.datum() > reg.zwoc()) {
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 ]);
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);
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 ]);
539 }
else if (reg.datum() < reg.zgoc()) {
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 ]);
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);
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 ]);
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);
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 ]);
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 ]);
589 template <
class Grid,
592 std::vector< std::vector<double> >
593 phasePressures(
const Grid& G,
595 const CellRange& cells,
598 std::array<double,2> span =
599 {{ std::numeric_limits<double>::max() ,
600 -std::numeric_limits<double>::max() }};
605 assert (UgGridHelpers::dimensions(G) == 3);
607 const int nd = UgGridHelpers::dimensions(G);
622 auto cell2Faces = UgGridHelpers::cell2Faces(G);
623 auto faceVertices = UgGridHelpers::face2Vertices(G);
625 for (
typename CellRange::const_iterator
626 ci = cells.begin(), ce = cells.end();
627 ci != ce; ++ci, ++ncell)
629 for (
auto fi=cell2Faces[*ci].begin(),
630 fe=cell2Faces[*ci].end();
634 for (
auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end();
637 const double z = UgGridHelpers::vertexCoordinates(G, *i)[nd-1];
639 if (z < span[0]) { span[0] = z; }
640 if (z > span[1]) { span[1] = z; }
645 const int np = reg.phaseUsage().num_phases;
647 typedef std::vector<double> pval;
648 std::vector<pval> press(np, pval(ncell, 0.0));
650 const double zwoc = reg.zwoc ();
651 const double zgoc = reg.zgoc ();
654 span[0] = std::min(span[0],zgoc);
655 span[1] = std::max(span[1],zwoc);
657 Details::equilibrateOWG(G, reg, grav, span, cells, press);
662 template <
class Grid,
666 temperature(
const Grid& ,
668 const CellRange& cells)
671 return std::vector<double>(cells.size(), 273.15 + 20.0);
674 template <
class Gr
id,
class Region,
class CellRange>
675 std::vector< std::vector<double> >
676 phaseSaturations(
const Grid& G,
678 const CellRange& cells,
679 BlackoilPropertiesInterface& props,
680 const std::vector<double> swat_init,
681 std::vector< std::vector<double> >& phase_pressures)
683 if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
684 OPM_THROW(std::runtime_error,
"Cannot initialise: not handling water-gas cases.");
687 std::vector< std::vector<double> > phase_saturations = phase_pressures;
688 double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
689 double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
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);
705 const double cellDepth = UgGridHelpers::cellCenterDepth(G,
707 sw =
satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,
false);
708 phase_saturations[waterpos][local_index] = sw;
711 const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
712 if (swat_init.empty()) {
713 sw =
satFromPc(props, waterpos, cell, pcov);
714 phase_saturations[waterpos][local_index] = sw;
716 sw = swat_init[cell];
717 props.swatInitScaling(cell, pcov, sw);
718 phase_saturations[waterpos][local_index] = sw;
725 const double cellDepth = UgGridHelpers::cellCenterDepth(G,
727 sg =
satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,
true);
728 phase_saturations[gaspos][local_index] = sg;
732 const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
733 const double increasing =
true;
734 sg =
satFromPc(props, gaspos, cell, pcog, increasing);
735 phase_saturations[gaspos][local_index] = sg;
738 if (gas && water && (sg + sw > 1.0)) {
743 const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
744 if (! swat_init.empty()) {
748 props.swatInitScaling(cell, pcgw, sw);
752 phase_saturations[waterpos][local_index] = sw;
753 phase_saturations[gaspos][local_index] = sg;
755 double pc[BlackoilPhases::MaxNumPhases];
756 double sat[BlackoilPhases::MaxNumPhases];
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];
763 phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
766 double pc[BlackoilPhases::MaxNumPhases];
767 double sat[BlackoilPhases::MaxNumPhases];
768 double threshold_sat = 1.0e-6;
772 sat[waterpos] = smax[waterpos];
773 sat[oilpos] -= sat[waterpos];
776 sat[gaspos] = smax[gaspos];
777 sat[oilpos] -= sat[gaspos];
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];
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];
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];
799 return phase_saturations;
821 template <
class Gr
id,
class CellRangeType>
823 const CellRangeType& cells,
824 const std::vector<double> oil_pressure,
825 const std::vector<double>& temperature,
827 const std::vector<double> gas_saturation)
829 assert(UgGridHelpers::dimensions(grid) == 3);
830 std::vector<double> rs(cells.size());
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]);
847 inline std::vector<double>
848 convertSats(
const std::vector< std::vector<double> >& sat)
850 const auto np = sat.size();
851 const auto nc = sat[0].size();
853 std::vector<double> s(np * nc);
855 for (decltype(sat.size()) p = 0; p < np; ++p) {
856 const auto& sat_p = sat[p];
857 double* sp = & s[0*nc + p];
859 for (decltype(sat[0].size()) c = 0;
860 c < nc; ++c, sp += np)
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)
898 typedef EQUIL::DeckDependent::InitialStateComputer ISC;
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];
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();
923 initBlackoilSurfvolUsingRSorRV(UgGridHelpers::numCells(grid), props, state);
930 #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
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: 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