22 #ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
23 #define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
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));
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>
116 const double norm_grav)
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>
158 const double norm_grav)
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>
208 const double norm_grav)
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 {
270 namespace PhaseIndex {
305 namespace PhasePressure {
306 template <
class Grid,
311 const std::array<PressFunction, 2>& f ,
313 const CellRange& cells,
314 std::vector<double>& p )
318 enum { up = 0, down = 1 };
320 std::vector<double>::size_type c = 0;
321 for (
typename CellRange::const_iterator
322 ci = cells.begin(), ce = cells.end();
325 assert (c < p.size());
328 p[c] = (z <
split) ? f[up](z) : f[down](z);
332 template <
class Grid,
338 const std::array<double,2>& span ,
341 const CellRange& cells ,
342 std::vector<double>& press )
345 typedef Water<typename Region::CalcDensity> ODE;
350 const double T = 273.15 + 20;
351 ODE drho(T, reg.densityCalculator(), pu.
num_phases, wix, grav);
355 if (reg.datum() > reg.zwoc()) {
360 p0 = po_woc - reg.pcow_woc();
363 std::array<double,2> up = {{ z0, span[0] }};
364 std::array<double,2> down = {{ z0, span[1] }};
367 std::array<WPress,2> wpress = {
369 WPress(drho, up , p0, 100)
371 WPress(drho, down, p0, 100)
375 assign(G, wpress, z0, cells, press);
377 if (reg.datum() > reg.zwoc()) {
379 po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc();
383 template <
class Grid,
389 const std::array<double,2>& span ,
391 const CellRange& cells ,
392 std::vector<double>& press ,
397 typedef Oil<
typename Region::CalcDensity,
398 typename Region::CalcDissolution> ODE;
404 const double T = 273.15 + 20;
406 reg.densityCalculator(),
407 reg.dissolutionCalculator(),
412 if (reg.datum() > reg.zwoc()) {
415 }
else if (reg.datum() < reg.zgoc()) {
423 std::array<double,2> up = {{ z0, span[0] }};
424 std::array<double,2> down = {{ z0, span[1] }};
427 std::array<OPress,2> opress = {
429 OPress(drho, up , p0, 100)
431 OPress(drho, down, p0, 100)
435 assign(G, opress, z0, cells, press);
437 const double woc = reg.zwoc();
438 if (z0 > woc) { po_woc = opress[0](woc); }
439 else if (z0 < woc) { po_woc = opress[1](woc); }
440 else { po_woc = p0; }
442 const double goc = reg.zgoc();
443 if (z0 > goc) { po_goc = opress[0](goc); }
444 else if (z0 < goc) { po_goc = opress[1](goc); }
445 else { po_goc = p0; }
448 template <
class Grid,
454 const std::array<double,2>& span ,
457 const CellRange& cells ,
458 std::vector<double>& press )
461 typedef Gas<
typename Region::CalcDensity,
462 typename Region::CalcEvaporation> ODE;
469 const double T = 273.15 + 20;
471 reg.densityCalculator(),
472 reg.evaporationCalculator(),
477 if (reg.datum() < reg.zgoc()) {
482 p0 = po_goc + reg.pcgo_goc();
485 std::array<double,2> up = {{ z0, span[0] }};
486 std::array<double,2> down = {{ z0, span[1] }};
489 std::array<GPress,2> gpress = {
491 GPress(drho, up , p0, 100)
493 GPress(drho, down, p0, 100)
497 assign(G, gpress, z0, cells, press);
499 if (reg.datum() < reg.zgoc()) {
501 po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc();
506 template <
class Grid,
513 const std::array<double,2>& span,
514 const CellRange& cells,
515 std::vector< std::vector<double> >& press)
519 if (reg.datum() > reg.zwoc()) {
526 cells, press[ wix ]);
532 press[ oix ], po_woc, po_goc);
538 cells, press[ gix ]);
540 }
else if (reg.datum() < reg.zgoc()) {
547 cells, press[ gix ]);
553 press[ oix ], po_woc, po_goc);
559 cells, press[ wix ]);
568 press[ oix ], po_woc, po_goc);
574 cells, press[ wix ]);
580 cells, press[ gix ]);
590 template <
class Grid,
593 std::vector< std::vector<double> >
596 const CellRange& cells,
599 std::array<double,2> span =
600 {{ std::numeric_limits<double>::max() ,
601 -std::numeric_limits<double>::max() }};
626 for (
typename CellRange::const_iterator
627 ci = cells.begin(), ce = cells.end();
628 ci != ce; ++ci, ++ncell)
635 for (
auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end();
640 if (z < span[0]) { span[0] = z; }
641 if (z > span[1]) { span[1] = z; }
646 const int np = reg.phaseUsage().num_phases;
648 typedef std::vector<double> pval;
649 std::vector<pval> press(np, pval(ncell, 0.0));
651 const double zwoc = reg.zwoc ();
652 const double zgoc = reg.zgoc ();
655 span[0] = std::min(span[0],zgoc);
656 span[1] = std::max(span[1],zwoc);
663 template <
class Grid,
669 const CellRange& cells)
672 return std::vector<double>(cells.size(), 273.15 + 20.0);
675 template <
class Gr
id,
class Region,
class CellRange>
676 std::vector< std::vector<double> >
679 const CellRange& cells,
681 const std::vector<double> swat_init,
682 std::vector< std::vector<double> >& phase_pressures)
685 OPM_THROW(std::runtime_error,
"Cannot initialise: not handling water-gas cases.");
688 std::vector< std::vector<double> > phase_saturations = phase_pressures;
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);
710 sw =
satFromDepth(props,cellDepth,reg.zwoc(),waterpos,cell,
false);
711 phase_saturations[waterpos][local_index] = sw;
714 const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
715 if (swat_init.empty()) {
716 sw =
satFromPc(props, waterpos, cell, pcov);
717 phase_saturations[waterpos][local_index] = sw;
719 sw = swat_init[cell];
721 phase_saturations[waterpos][local_index] = sw;
732 sg =
satFromDepth(props,cellDepth,reg.zgoc(),gaspos,cell,
true);
733 phase_saturations[gaspos][local_index] = sg;
737 const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
738 const double increasing =
true;
739 sg =
satFromPc(props, gaspos, cell, pcog, increasing);
740 phase_saturations[gaspos][local_index] = sg;
743 if (gas && water && (sg + sw > 1.0)) {
748 const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
749 if (! swat_init.empty()) {
757 phase_saturations[waterpos][local_index] = sw;
758 phase_saturations[gaspos][local_index] = 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];
768 phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
773 double threshold_sat = 1.0e-6;
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];
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];
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];
798 return phase_saturations;
820 template <
class Gr
id,
class CellRangeType>
822 const CellRangeType& cells,
823 const std::vector<double> oil_pressure,
826 const std::vector<double> gas_saturation)
829 std::vector<double> rs(cells.size());
831 for (
auto it = cells.begin(); it != cells.end(); ++it, ++count) {
833 rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]);
846 inline std::vector<double>
849 const auto np = sat.size();
850 const auto nc = sat[0].size();
852 std::vector<double> s(np * nc);
854 for (decltype(sat.size()) p = 0; p < np; ++p) {
855 const auto& sat_p = sat[p];
856 double* sp = & s[0*nc + p];
858 for (decltype(sat[0].size()) c = 0;
859 c < nc; ++c, sp += np)
888 const Opm::DeckConstPtr deck,
889 const Opm::EclipseStateConstPtr eclipseState,
893 typedef Equil::DeckDependent::InitialStateComputer ISC;
894 ISC isc(props, deck, eclipseState, grid, gravity);
899 state.
pressure() = isc.press()[ref_phase];
902 state.
rv() = isc.rv();
910 #endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
Definition: initStateEquil_impl.hpp:110
int oil(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:283
int gas(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:294
std::pair< std::string, std::string > split(const std::string &name)
Water(const double temp, const Density &rho, const int np, const int ix, const double norm_grav)
Definition: initStateEquil_impl.hpp:112
void initStateEquil(const Grid &grid, const BlackoilPropertiesInterface &props, const Opm::DeckConstPtr deck, const Opm::EclipseStateConstPtr eclipseState, const double gravity, BlackoilState &state)
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:452
double operator()(const double depth, const double press) const
Definition: initStateEquil_impl.hpp:171
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:510
Definition: AnisotropicEikonal.hpp:43
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
Definition: BlackoilPhases.hpp:32
std::vector< double > & rv()
Definition: BlackoilState.hpp:55
bool gas(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:264
std::vector< double > & saturation()
Definition: SimulatorState.hpp:60
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
Definition: initStateEquil_impl.hpp:150
Face2VerticesTraits< UnstructuredGrid >::Type face2Vertices(const UnstructuredGrid &grid)
Get the face to vertices mapping of a grid.
virtual void capPress(const int n, const double *s, const int *cells, double *pc, double *dpcds) const =0
int phase_used[MaxNumPhases]
Definition: BlackoilPhases.hpp:39
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
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:336
bool isConstPc(const BlackoilPropertiesInterface &props, const int phase, const int cell)
Return true if capillary pressure function is constant.
Definition: EquilibrationHelpers.hpp:890
bool water(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:252
int numCells(const UnstructuredGrid &grid)
Get the number of cells of a grid.
Definition: initStateEquil_impl.hpp:200
Definition: BlackoilPropertiesInterface.hpp:37
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:387
int phase_pos[MaxNumPhases]
Definition: BlackoilPhases.hpp:40
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
std::vector< std::vector< double > > phasePressures(const Grid &G, const Region ®, const CellRange &cells, const double grav)
Definition: initStateEquil_impl.hpp:594
virtual void satRange(const int n, const int *cells, double *smin, double *smax) const =0
int num_phases
Definition: BlackoilPhases.hpp:38
Simulator state for a blackoil simulator.
Definition: BlackoilState.hpp:33
double cellCentroidCoordinate(const UnstructuredGrid &grid, int cell_index, int coordinate)
Get a coordinate of a specific cell centroid.
int dimensions(const UnstructuredGrid &grid)
Get the dimensions of a grid.
std::vector< double > convertSats(const std::vector< std::vector< double > > &sat)
Definition: initStateEquil_impl.hpp:847
virtual PhaseUsage phaseUsage() const =0
double satFromPc(const BlackoilPropertiesInterface &props, const int phase, const int cell, const double target_pc, const bool increasing=false)
Definition: EquilibrationHelpers.hpp:760
const double gravity
Definition: Units.hpp:120
double operator()(const double x) const
Definition: initStateEquil_impl.hpp:75
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
Definition: BlackoilPhases.hpp:32
Definition: initStateEquil_impl.hpp:39
double operator()(const double depth, const double press) const
Definition: initStateEquil_impl.hpp:221
const UnstructuredGrid & grid
Definition: ColumnExtract.hpp:31
static const int MaxNumPhases
Definition: BlackoilPhases.hpp:30
bool oil(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:258
Cell2FacesTraits< UnstructuredGrid >::Type cell2Faces(const UnstructuredGrid &grid)
Get the cell to faces mapping of a grid.
std::vector< double > & gasoilratio()
Definition: BlackoilState.hpp:54
void initBlackoilSurfvolUsingRSorRV(const UnstructuredGrid &grid, const Props &props, State &state)
Definition: initState_impl.hpp:724
int water(const PhaseUsage &pu)
Definition: initStateEquil_impl.hpp:272
virtual void swatInitScaling(const int cell, const double pcow, double &swat)=0
RK4IVP(const RHS &f, const std::array< double, 2 > &span, const double y0, const int N)
Definition: initStateEquil_impl.hpp:41
double operator()(const double, const double press) const
Definition: initStateEquil_impl.hpp:127
std::vector< double > & pressure()
Definition: SimulatorState.hpp:56
const double * vertexCoordinates(const UnstructuredGrid &grid, int index)
Get the coordinates of a vertex of the grid.
Definition: BlackoilPhases.hpp:32
std::vector< double > temperature(const Grid &, const Region &, const CellRange &cells)
Definition: initStateEquil_impl.hpp:667
Definition: BlackoilPhases.hpp:36
Definition: EquilibrationHelpers.hpp:165
double satFromSumOfPcs(const BlackoilPropertiesInterface &props, const int phase1, const int phase2, const int cell, const double target_pc)
Definition: EquilibrationHelpers.hpp:835
std::vector< std::vector< double > > phaseSaturations(const Grid &G, const Region ®, const CellRange &cells, BlackoilPropertiesInterface &props, const std::vector< double > swat_init, std::vector< std::vector< double > > &phase_pressures)
Definition: initStateEquil_impl.hpp:677