diff --git a/ebos/equil/initStateEquil.hpp b/ebos/equil/initStateEquil.hpp index ae7fc0efd..e77f25ab3 100644 --- a/ebos/equil/initStateEquil.hpp +++ b/ebos/equil/initStateEquil.hpp @@ -56,391 +56,1079 @@ #include #include -namespace Ewoms +namespace Ewoms { +/** + * Types and routines that collectively implement a basic + * ECLIPSE-style equilibration-based initialisation scheme. + * + * This namespace is intentionally nested to avoid name clashes + * with other parts of OPM. + */ +namespace EQUIL { +namespace Details { +template +class RK4IVP { +public: + RK4IVP(const RHS& f , + const std::array& span, + const double y0 , + const int N ) + : N_(N) + , span_(span) + { + const double h = stepsize(); + const double h2 = h / 2; + const double h6 = h / 6; + + y_.reserve(N + 1); + f_.reserve(N + 1); + + y_.push_back(y0); + f_.push_back(f(span_[0], y0)); + + for (int i = 0; i < N; ++i) { + const double x = span_[0] + i*h; + const double y = y_.back(); + + const double k1 = f_[i]; + const double k2 = f(x + h2, y + h2*k1); + const double k3 = f(x + h2, y + h2*k2); + const double k4 = f(x + h , y + h*k3); + + y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); + f_.push_back(f(x + h, y_.back())); + } + + assert (y_.size() == std::vector::size_type(N + 1)); + } + + double + operator()(const double x) const + { + // Dense output (O(h**3)) according to Shampine + // (Hermite interpolation) + const double h = stepsize(); + int i = (x - span_[0]) / h; + const double t = (x - (span_[0] + i*h)) / h; + + // Crude handling of evaluation point outside "span_"; + if (i < 0) { i = 0; } + if (N_ <= i) { i = N_ - 1; } + + const double y0 = y_[i], y1 = y_[i + 1]; + const double f0 = f_[i], f1 = f_[i + 1]; + + double u = (1 - 2*t) * (y1 - y0); + u += h * ((t - 1)*f0 + t*f1); + u *= t * (t - 1); + u += (1 - t)*y0 + t*y1; + + return u; + } + +private: + int N_; + std::array span_; + std::vector y_; + std::vector f_; + + double + stepsize() const { return (span_[1] - span_[0]) / N_; } +}; + +namespace PhasePressODE { +template +class Water { +public: + Water(const double temp, + const int pvtRegionIdx, + const double norm_grav) + : temp_(temp) + , pvtRegionIdx_(pvtRegionIdx) + , g_(norm_grav) + { + } + + double + operator()(const double /* depth */, + const double press) const + { + return this->density(press) * g_; + } + +private: + const double temp_; + const int pvtRegionIdx_; + const double g_; + + double + density(const double press) const + { + double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); + return rho; + } +}; + +template +class Oil { +public: + Oil(const double temp, + const RS& rs, + const int pvtRegionIdx, + const double norm_grav) + : temp_(temp) + , rs_(rs) + , pvtRegionIdx_(pvtRegionIdx) + , g_(norm_grav) + { + } + + double + operator()(const double depth, + const double press) const + { + return this->density(depth, press) * g_; + } + +private: + const double temp_; + const RS& rs_; + const int pvtRegionIdx_; + const double g_; + + double + density(const double depth, + const double press) const + { + double rs = rs_(depth, press, temp_); + double bOil = 0.0; + if ( !FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press) ) { + bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + } else { + bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs); + } + double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + if (FluidSystem::enableDissolvedGas()) { + rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + } + + return rho; + } +}; + +template +class Gas { +public: + Gas(const double temp, + const RV& rv, + const int pvtRegionIdx, + const double norm_grav) + : temp_(temp) + , rv_(rv) + , pvtRegionIdx_(pvtRegionIdx) + , g_(norm_grav) + { + } + + double + operator()(const double depth, + const double press) const + { + return this->density(depth, press) * g_; + } + +private: + const double temp_; + const RV& rv_; + const int pvtRegionIdx_; + const double g_; + + double + density(const double depth, + const double press) const + { + double rv = rv_(depth, press, temp_); + double bGas = 0.0; + if ( !FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) ) { + bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); + } else { + bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv); + } + double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); + if (FluidSystem::enableVaporizedOil()) { + rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); + } + + return rho; + } +}; +} // namespace PhasePressODE + + +namespace PhasePressure { +template +void +assign(const Grid& grid , + const std::array& f , + const double split, + const CellRange& cells, + std::vector& p ) { + enum { up = 0, down = 1 }; - /** - * Types and routines that collectively implement a basic - * ECLIPSE-style equilibration-based initialisation scheme. - * - * This namespace is intentionally nested to avoid name clashes - * with other parts of OPM. - */ - namespace EQUIL { + std::vector::size_type c = 0; + for (typename CellRange::const_iterator + ci = cells.begin(), ce = cells.end(); + ci != ce; ++ci, ++c) + { + assert (c < p.size()); - /** - * Compute initial phase pressures by means of equilibration. - * - * This function uses the information contained in an - * equilibration record (i.e., depths and pressurs) as well as - * a density calculator and related data to vertically - * integrate the phase pressure ODE - * \f[ - * \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} = - * \rho_{\alpha}(z,p_{\alpha})\cdot g - * \f] - * in which \f$\rho_{\alpha}$ denotes the fluid density of - * fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the - * corresponding phase pressure, \f$z\f$ is the depth and - * \f$g\f$ is the acceleration due to gravity (assumed - * directed downwords, in the positive \f$z\f$ direction). - * - * \tparam Region Type of an equilibration region information - * base. Typically an instance of the EquilReg - * class template. - * - * \tparam CellRange Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \param[in] G Grid. - * \param[in] reg Current equilibration region. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] grav Acceleration of gravity. - * - * \return Phase pressures, one vector for each active phase, - * of pressure values in each cell in the current - * equilibration region. - */ - template - std::vector< std::vector > - phasePressures(const Grid& G, - const Region& reg, - const CellRange& cells, - const double grav = Opm::unit::gravity); + const double z = Opm::UgGridHelpers::cellCenterDepth(grid, *ci); + p[c] = (z < split) ? f[up](z) : f[down](z); + } +} +template < class FluidSystem, + class Grid, + class Region, + class CellRange> +void +water(const Grid& grid , + const Region& reg , + const std::array& span , + const double grav , + double& po_woc, + const CellRange& cells , + std::vector& press ) +{ + using PhasePressODE::Water; + typedef Water ODE; + const double T = 273.15 + 20; // standard temperature for now + ODE drho(T, reg.pvtIdx() , grav); - /** - * Compute initial phase saturations by means of equilibration. - * - * \tparam FluidSystem The FluidSystem from opm-material - * Must be initialized before used. - * - * \tparam Grid Type of the grid - * - * \tparam Region Type of an equilibration region information - * base. Typically an instance of the EquilReg - * class template. - * - * \tparam CellRange Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \tparam MaterialLawManager The MaterialLawManager from opm-material - * - * \param[in] G Grid. - * \param[in] reg Current equilibration region. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] materialLawManager The MaterialLawManager from opm-material - * \param[in] swat_init A vector of initial water saturations. - * The capillary pressure is scaled to fit these values - * \param[in] phase_pressures Phase pressures, one vector for each active phase, - * of pressure values in each cell in the current - * equilibration region. - * \return Phase saturations, one vector for each phase, each containing - * one saturation value per cell in the region. - */ - template - std::vector< std::vector > - phaseSaturations(const Grid& grid, - const Region& reg, - const CellRange& cells, - MaterialLawManager& materialLawManager, - const std::vector swat_init, - std::vector< std::vector >& phase_pressures); + double z0; + double p0; + if (reg.datum() > reg.zwoc()) {//Datum in water zone + z0 = reg.datum(); + p0 = reg.pressure(); + } else { + z0 = reg.zwoc(); + p0 = po_woc - reg.pcow_woc(); // Water pressure at contact + } + std::array up = {{ z0, span[0] }}; + std::array down = {{ z0, span[1] }}; + typedef Details::RK4IVP WPress; + std::array wpress = { + { + WPress(drho, up , p0, 2000) + , + WPress(drho, down, p0, 2000) + } + }; - /** - * Compute initial Rs values. - * - * \tparam CellRangeType Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \param[in] grid Grid. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] oil_pressure Oil pressure for each cell in range. - * \param[in] temperature Temperature for each cell in range. - * \param[in] rs_func Rs as function of pressure and depth. - * \return Rs values, one for each cell in the 'cells' range. - */ - template - std::vector computeRs(const Grid& grid, - const CellRangeType& cells, - const std::vector oil_pressure, - const std::vector& temperature, - const Miscibility::RsFunction& rs_func, - const std::vector gas_saturation); + assign(grid, wpress, z0, cells, press); - namespace DeckDependent { - inline - std::vector - getEquil(const Opm::EclipseState& state) + if (reg.datum() > reg.zwoc()) { + // Return oil pressure at contact + po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc(); + } +} + +template +void +oil(const Grid& grid , + const Region& reg , + const std::array& span , + const double grav , + const CellRange& cells , + std::vector& press , + double& po_woc, + double& po_goc) +{ + using PhasePressODE::Oil; + typedef Oil ODE; + + const double T = 273.15 + 20; // standard temperature for now + ODE drho(T, reg.dissolutionCalculator(), + reg.pvtIdx(), grav); + + double z0; + double p0; + if (reg.datum() > reg.zwoc()) {//Datum in water zone, po_woc given + z0 = reg.zwoc(); + p0 = po_woc; + } else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, po_goc given + z0 = reg.zgoc(); + p0 = po_goc; + } else { //Datum in oil zone + z0 = reg.datum(); + p0 = reg.pressure(); + } + + std::array up = {{ z0, span[0] }}; + std::array down = {{ z0, span[1] }}; + + typedef Details::RK4IVP OPress; + std::array opress = { + { + OPress(drho, up , p0, 2000) + , + OPress(drho, down, p0, 2000) + } + }; + + assign(grid, opress, z0, cells, press); + + const double woc = reg.zwoc(); + if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum + else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum + else { po_woc = p0; } // WOC *at* datum + + const double goc = reg.zgoc(); + if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum + else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum + else { po_goc = p0; } // GOC *at* datum +} + +template +void +gas(const Grid& grid , + const Region& reg , + const std::array& span , + const double grav , + double& po_goc, + const CellRange& cells , + std::vector& press ) +{ + using PhasePressODE::Gas; + typedef Gas ODE; + + const double T = 273.15 + 20; // standard temperature for now + ODE drho(T, reg.evaporationCalculator(), + reg.pvtIdx(), grav); + + double z0; + double p0; + if (reg.datum() < reg.zgoc()) {//Datum in gas zone + z0 = reg.datum(); + p0 = reg.pressure(); + } else { + z0 = reg.zgoc(); + p0 = po_goc + reg.pcgo_goc(); // Gas pressure at contact + } + + std::array up = {{ z0, span[0] }}; + std::array down = {{ z0, span[1] }}; + + typedef Details::RK4IVP GPress; + std::array gpress = { + { + GPress(drho, up , p0, 2000) + , + GPress(drho, down, p0, 2000) + } + }; + + assign(grid, gpress, z0, cells, press); + + if (reg.datum() < reg.zgoc()) { + // Return oil pressure at contact + po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc(); + } +} +} // namespace PhasePressure + +template +void +equilibrateOWG(const Grid& grid, + const Region& reg, + const double grav, + const std::array& span, + const CellRange& cells, + std::vector< std::vector >& press) +{ + const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); + const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + const int oilpos = FluidSystem::oilPhaseIdx; + const int waterpos = FluidSystem::waterPhaseIdx; + const int gaspos = FluidSystem::gasPhaseIdx; + + if (reg.datum() > reg.zwoc()) { // Datum in water zone + double po_woc = -1; + double po_goc = -1; + + if (water) { + PhasePressure::water(grid, reg, span, grav, po_woc, + cells, press[ waterpos ]); + } + + if (oil) { + PhasePressure::oil(grid, reg, span, grav, cells, + press[ oilpos ], po_woc, po_goc); + } + + if (gas) { + PhasePressure::gas(grid, reg, span, grav, po_goc, + cells, press[ gaspos ]); + } + } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone + double po_woc = -1; + double po_goc = -1; + + if (gas) { + PhasePressure::gas(grid, reg, span, grav, po_goc, + cells, press[ gaspos ]); + } + + if (oil) { + PhasePressure::oil(grid, reg, span, grav, cells, + press[ oilpos ], po_woc, po_goc); + } + + if (water) { + PhasePressure::water(grid, reg, span, grav, po_woc, + cells, press[ waterpos ]); + } + } else { // Datum in oil zone + double po_woc = -1; + double po_goc = -1; + + if (oil) { + PhasePressure::oil(grid, reg, span, grav, cells, + press[ oilpos ], po_woc, po_goc); + } + + if (water) { + PhasePressure::water(grid, reg, span, grav, po_woc, + cells, press[ waterpos ]); + } + + if (gas) { + PhasePressure::gas(grid, reg, span, grav, po_goc, + cells, press[ gaspos ]); + } + } +} +} // namespace Details + +/** + * Compute initial phase pressures by means of equilibration. + * + * This function uses the information contained in an + * equilibration record (i.e., depths and pressurs) as well as + * a density calculator and related data to vertically + * integrate the phase pressure ODE + * \f[ + * \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} = + * \rho_{\alpha}(z,p_{\alpha})\cdot g + * \f] + * in which \f$\rho_{\alpha}$ denotes the fluid density of + * fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the + * corresponding phase pressure, \f$z\f$ is the depth and + * \f$g\f$ is the acceleration due to gravity (assumed + * directed downwords, in the positive \f$z\f$ direction). + * + * \tparam Region Type of an equilibration region information + * base. Typically an instance of the EquilReg + * class template. + * + * \tparam CellRange Type of cell range that demarcates the + * cells pertaining to the current + * equilibration region. Must implement + * methods begin() and end() to bound the range + * as well as provide an inner type, + * const_iterator, to traverse the range. + * + * \param[in] grid Grid. + * \param[in] reg Current equilibration region. + * \param[in] cells Range that spans the cells of the current + * equilibration region. + * \param[in] grav Acceleration of gravity. + * + * \return Phase pressures, one vector for each active phase, + * of pressure values in each cell in the current + * equilibration region. + */ +template +std::vector< std::vector > +phasePressures(const Grid& grid, + const Region& reg, + const CellRange& cells, + const double grav = Opm::unit::gravity) +{ + std::array span = + {{ std::numeric_limits::max() , + -std::numeric_limits::max() }}; // Symm. about 0. + + int ncell = 0; + { + // This code is only supported in three space dimensions + assert (Opm::UgGridHelpers::dimensions(grid) == 3); + + const int nd = Opm::UgGridHelpers::dimensions(grid); + + // Define vertical span as + // + // [minimum(node depth(cells)), maximum(node depth(cells))] + // + // Note: We use a sledgehammer approach--looping all + // the nodes of all the faces of all the 'cells'--to + // compute those bounds. This necessarily entails + // visiting some nodes (and faces) multiple times. + // + // Note: The implementation of 'RK4IVP<>' implicitly + // imposes the requirement that cell centroids are all + // within this vertical span. That requirement is not + // checked. + auto cell2Faces = Opm::UgGridHelpers::cell2Faces(grid); + auto faceVertices = Opm::UgGridHelpers::face2Vertices(grid); + + for (typename CellRange::const_iterator + ci = cells.begin(), ce = cells.end(); + ci != ce; ++ci, ++ncell) + { + for (auto fi=cell2Faces[*ci].begin(), + fe=cell2Faces[*ci].end(); + fi != fe; + ++fi) { - const auto& init = state.getInitConfig(); - - if( !init.hasEquil() ) { - OPM_THROW(std::domain_error, "Deck does not provide equilibration data."); - } - - const auto& equil = init.getEquil(); - return { equil.begin(), equil.end() }; - } - - template - inline - std::vector - equilnum(const Opm::EclipseState& eclipseState, - const Grid& G ) - { - std::vector eqlnum; - if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) { - const int nc = Opm::UgGridHelpers::numCells(G); - eqlnum.resize(nc); - const std::vector& e = - eclipseState.get3DProperties().getIntGridProperty("EQLNUM").getData(); - const int* gc = Opm::UgGridHelpers::globalCell(G); - for (int cell = 0; cell < nc; ++cell) { - const int deck_pos = (gc == NULL) ? cell : gc[cell]; - eqlnum[cell] = e[deck_pos] - 1; - } - } - else { - // No explicit equilibration region. - // All cells in region zero. - eqlnum.assign(Opm::UgGridHelpers::numCells(G), 0); - } - - return eqlnum; - } - - template - class InitialStateComputer { - public: - template - InitialStateComputer(MaterialLawManager& materialLawManager, - const Opm::EclipseState& eclipseState, - const Grid& G , - const double grav = Opm::unit::gravity, - const bool applySwatInit = true - ) - : pp_(FluidSystem::numPhases, - std::vector(Opm::UgGridHelpers::numCells(G))), - sat_(FluidSystem::numPhases, - std::vector(Opm::UgGridHelpers::numCells(G))), - rs_(Opm::UgGridHelpers::numCells(G)), - rv_(Opm::UgGridHelpers::numCells(G)) + for (auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end(); + i != e; ++i) { - //Check for presence of kw SWATINIT - if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) { - const std::vector& swat_init_ecl = eclipseState. - get3DProperties().getDoubleGridProperty("SWATINIT").getData(); - const int nc = Opm::UgGridHelpers::numCells(G); - swat_init_.resize(nc); - const int* gc = Opm::UgGridHelpers::globalCell(G); - for (int c = 0; c < nc; ++c) { - const int deck_pos = (gc == NULL) ? c : gc[c]; - swat_init_[c] = swat_init_ecl[deck_pos]; - } + const double z = Opm::UgGridHelpers::vertexCoordinates(grid, *i)[nd-1]; + + if (z < span[0]) { span[0] = z; } + if (z > span[1]) { span[1] = z; } + } + } + } + } + const int np = FluidSystem::numPhases; //reg.phaseUsage().num_phases; + + typedef std::vector pval; + std::vector press(np, pval(ncell, 0.0)); + + const double zwoc = reg.zwoc (); + const double zgoc = reg.zgoc (); + + // make sure goc and woc is within the span for the phase pressure calculation + span[0] = std::min(span[0],zgoc); + span[1] = std::max(span[1],zwoc); + + Details::equilibrateOWG(grid, reg, grav, span, cells, press); + + return press; +} + +template +std::vector +temperature(const Grid& /* G */, + const Region& /* reg */, + const CellRange& cells) +{ + // use the standard temperature for everything for now + return std::vector(cells.size(), 273.15 + 20.0); +} + +/** + * Compute initial phase saturations by means of equilibration. + * + * \tparam FluidSystem The FluidSystem from opm-material + * Must be initialized before used. + * + * \tparam Grid Type of the grid + * + * \tparam Region Type of an equilibration region information + * base. Typically an instance of the EquilReg + * class template. + * + * \tparam CellRange Type of cell range that demarcates the + * cells pertaining to the current + * equilibration region. Must implement + * methods begin() and end() to bound the range + * as well as provide an inner type, + * const_iterator, to traverse the range. + * + * \tparam MaterialLawManager The MaterialLawManager from opm-material + * + * \param[in] grid Grid. + * \param[in] reg Current equilibration region. + * \param[in] cells Range that spans the cells of the current + * equilibration region. + * \param[in] materialLawManager The MaterialLawManager from opm-material + * \param[in] swat_init A vector of initial water saturations. + * The capillary pressure is scaled to fit these values + * \param[in] phase_pressures Phase pressures, one vector for each active phase, + * of pressure values in each cell in the current + * equilibration region. + * \return Phase saturations, one vector for each phase, each containing + * one saturation value per cell in the region. + */ +template +std::vector< std::vector > +phaseSaturations(const Grid& grid, + const Region& reg, + const CellRange& cells, + MaterialLawManager& materialLawManager, + const std::vector swat_init, + std::vector< std::vector >& phase_pressures) +{ + if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { + OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases."); + } + + std::vector< std::vector > phase_saturations = phase_pressures; // Just to get the right size. + + // Adjust oil pressure according to gas saturation and cap pressure + typedef Opm::SimpleModularFluidState SatOnlyFluidState; + + SatOnlyFluidState fluidState; + typedef typename MaterialLawManager::MaterialLaw MaterialLaw; + + const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); + const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + const int oilpos = FluidSystem::oilPhaseIdx; + const int waterpos = FluidSystem::waterPhaseIdx; + const int gaspos = FluidSystem::gasPhaseIdx; + std::vector::size_type local_index = 0; + for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) { + const int cell = *ci; + const auto& scaledDrainageInfo = + materialLawManager.oilWaterScaledEpsInfoDrainage(cell); + const auto& matParams = materialLawManager.materialLawParams(cell); + + // Find saturations from pressure differences by + // inverting capillary pressure functions. + double sw = 0.0; + if (water) { + if (isConstPc(materialLawManager,FluidSystem::waterPhaseIdx, cell)){ + const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(grid, + cell); + sw = satFromDepth(materialLawManager,cellDepth,reg.zwoc(),waterpos,cell,false); + phase_saturations[waterpos][local_index] = sw; + } + else{ + const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; + if (swat_init.empty()) { // Invert Pc to find sw + sw = satFromPc(materialLawManager, waterpos, cell, pcov); + phase_saturations[waterpos][local_index] = sw; + } else { // Scale Pc to reflect imposed sw + sw = swat_init[cell]; + sw = materialLawManager.applySwatinit(cell, pcov, sw); + phase_saturations[waterpos][local_index] = sw; + } + } + } + double sg = 0.0; + if (gas) { + if (isConstPc(materialLawManager,FluidSystem::gasPhaseIdx,cell)){ + const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(grid, + cell); + sg = satFromDepth(materialLawManager,cellDepth,reg.zgoc(),gaspos,cell,true); + phase_saturations[gaspos][local_index] = sg; + } + else{ + // Note that pcog is defined to be (pg - po), not (po - pg). + const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index]; + const double increasing = true; // pcog(sg) expected to be increasing function + sg = satFromPc(materialLawManager, gaspos, cell, pcog, increasing); + phase_saturations[gaspos][local_index] = sg; + } + } + if (gas && water && (sg + sw > 1.0)) { + // Overlapping gas-oil and oil-water transition + // zones can lead to unphysical saturations when + // treated as above. Must recalculate using gas-water + // capillary pressure. + const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index]; + if (! swat_init.empty()) { + // Re-scale Pc to reflect imposed sw for vanishing oil phase. + // This seems consistent with ecl, and fails to honour + // swat_init in case of non-trivial gas-oil cap pressure. + sw = materialLawManager.applySwatinit(cell, pcgw, sw); + } + sw = satFromSumOfPcs(materialLawManager, waterpos, gaspos, cell, pcgw); + sg = 1.0 - sw; + phase_saturations[waterpos][local_index] = sw; + phase_saturations[gaspos][local_index] = sg; + if ( water ) { + fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); + } + else { + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); + } + fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - sw - sg); + fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg); + + double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; + phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas; + } + phase_saturations[oilpos][local_index] = 1.0 - sw - sg; + + // Adjust phase pressures for max and min saturation ... + double threshold_sat = 1.0e-6; + + double so = 1.0; + double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 }; + if (water) { + double swu = scaledDrainageInfo.Swu; + fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu); + so -= swu; + } + if (gas) { + double sgu = scaledDrainageInfo.Sgu; + fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu); + so-= sgu; + } + fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); + + if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) { + fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu); + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; + phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat; + } else if (gas && sg > scaledDrainageInfo.Sgu-threshold_sat) { + fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu); + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; + phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas; + } + if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) { + fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl); + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; + phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pcGas; + } + if (water && sw < scaledDrainageInfo.Swl+threshold_sat) { + fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl); + MaterialLaw::capillaryPressures(pC, matParams, fluidState); + double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; + phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat; + } + } + return phase_saturations; +} + +/** + * Compute initial Rs values. + * + * \tparam CellRangeType Type of cell range that demarcates the + * cells pertaining to the current + * equilibration region. Must implement + * methods begin() and end() to bound the range + * as well as provide an inner type, + * const_iterator, to traverse the range. + * + * \param[in] grid Grid. + * \param[in] cells Range that spans the cells of the current + * equilibration region. + * \param[in] oil_pressure Oil pressure for each cell in range. + * \param[in] temperature Temperature for each cell in range. + * \param[in] rs_func Rs as function of pressure and depth. + * \return Rs values, one for each cell in the 'cells' range. + */ +template +std::vector computeRs(const Grid& grid, + const CellRangeType& cells, + const std::vector oil_pressure, + const std::vector& temperature, + const Miscibility::RsFunction& rs_func, + const std::vector gas_saturation) +{ + assert(Opm::UgGridHelpers::dimensions(grid) == 3); + std::vector rs(cells.size()); + int count = 0; + for (auto it = cells.begin(); it != cells.end(); ++it, ++count) { + const double depth = Opm::UgGridHelpers::cellCenterDepth(grid, *it); + rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]); + } + return rs; +} + +namespace DeckDependent { +inline +std::vector +getEquil(const Opm::EclipseState& state) +{ + const auto& init = state.getInitConfig(); + + if( !init.hasEquil() ) { + OPM_THROW(std::domain_error, "Deck does not provide equilibration data."); + } + + const auto& equil = init.getEquil(); + return { equil.begin(), equil.end() }; +} + +template +inline +std::vector +equilnum(const Opm::EclipseState& eclipseState, + const Grid& grid ) +{ + std::vector eqlnum; + if (eclipseState.get3DProperties().hasDeckIntGridProperty("EQLNUM")) { + const int nc = Opm::UgGridHelpers::numCells(grid); + eqlnum.resize(nc); + const std::vector& e = + eclipseState.get3DProperties().getIntGridProperty("EQLNUM").getData(); + const int* gc = Opm::UgGridHelpers::globalCell(grid); + for (int cell = 0; cell < nc; ++cell) { + const int deck_pos = (gc == NULL) ? cell : gc[cell]; + eqlnum[cell] = e[deck_pos] - 1; + } + } + else { + // No explicit equilibration region. + // All cells in region zero. + eqlnum.assign(Opm::UgGridHelpers::numCells(grid), 0); + } + + return eqlnum; +} + +template +class InitialStateComputer { +public: + template + InitialStateComputer(MaterialLawManager& materialLawManager, + const Opm::EclipseState& eclipseState, + const Grid& grid , + const double grav = Opm::unit::gravity, + const bool applySwatInit = true + ) + : pp_(FluidSystem::numPhases, + std::vector(Opm::UgGridHelpers::numCells(grid))), + sat_(FluidSystem::numPhases, + std::vector(Opm::UgGridHelpers::numCells(grid))), + rs_(Opm::UgGridHelpers::numCells(grid)), + rv_(Opm::UgGridHelpers::numCells(grid)) + { + //Check for presence of kw SWATINIT + if (eclipseState.get3DProperties().hasDeckDoubleGridProperty("SWATINIT") && applySwatInit) { + const std::vector& swat_init_ecl = eclipseState. + get3DProperties().getDoubleGridProperty("SWATINIT").getData(); + const int nc = Opm::UgGridHelpers::numCells(grid); + swat_init_.resize(nc); + const int* gc = Opm::UgGridHelpers::globalCell(grid); + for (int c = 0; c < nc; ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + swat_init_[c] = swat_init_ecl[deck_pos]; + } + } + // Get the equilibration records. + const std::vector rec = getEquil(eclipseState); + const auto& tables = eclipseState.getTableManager(); + // Create (inverse) region mapping. + const Ewoms::RegionMapping<> eqlmap(equilnum(eclipseState, grid)); + const int invalidRegion = -1; + regionPvtIdx_.resize(rec.size(), invalidRegion); + setRegionPvtIdx(grid, eclipseState, eqlmap); + + // Create Rs functions. + rs_func_.reserve(rec.size()); + if (FluidSystem::enableDissolvedGas()) { + const Opm::TableContainer& rsvdTables = tables.getRsvdTables(); + for (size_t i = 0; i < rec.size(); ++i) { + if (eqlmap.cells(i).empty()) + { + rs_func_.push_back(std::shared_ptr>()); + continue; + } + const int pvtIdx = regionPvtIdx_[i]; + if (!rec[i].liveOilInitConstantRs()) { + if (rsvdTables.size() <= 0 ) { + OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available."); } - // Get the equilibration records. - const std::vector rec = getEquil(eclipseState); - const auto& tables = eclipseState.getTableManager(); - // Create (inverse) region mapping. - const Ewoms::RegionMapping<> eqlmap(equilnum(eclipseState, grid)); - const int invalidRegion = -1; - regionPvtIdx_.resize(rec.size(), invalidRegion); - setRegionPvtIdx(G, eclipseState, eqlmap); - - // Create Rs functions. - rs_func_.reserve(rec.size()); - if (FluidSystem::enableDissolvedGas()) { - const Opm::TableContainer& rsvdTables = tables.getRsvdTables(); - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) - { - rs_func_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].liveOilInitConstantRs()) { - if (rsvdTables.size() <= 0 ) { - OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table not available."); - } - const Opm::RsvdTable& rsvdTable = rsvdTables.getTable(i); - std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); - rs_func_.push_back(std::make_shared>(pvtIdx, + const Opm::RsvdTable& rsvdTable = rsvdTables.getTable(i); + std::vector depthColumn = rsvdTable.getColumn("DEPTH").vectorCopy(); + std::vector rsColumn = rsvdTable.getColumn("RS").vectorCopy(); + rs_func_.push_back(std::make_shared>(pvtIdx, depthColumn , rsColumn)); - } else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - OPM_THROW(std::runtime_error, - "Cannot initialise: when no explicit RSVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); - } - const double p_contact = rec[i].datumDepthPressure(); - const double T_contact = 273.15 + 20; // standard temperature for now - rs_func_.push_back(std::make_shared>(pvtIdx, p_contact, T_contact)); - } - } - } else { - for (size_t i = 0; i < rec.size(); ++i) { - rs_func_.push_back(std::make_shared()); - } - } + } else { + if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { + OPM_THROW(std::runtime_error, + "Cannot initialise: when no explicit RSVD table is given, \n" + "datum depth must be at the gas-oil-contact. " + "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); + } + const double p_contact = rec[i].datumDepthPressure(); + const double T_contact = 273.15 + 20; // standard temperature for now + rs_func_.push_back(std::make_shared>(pvtIdx, p_contact, T_contact)); + } + } + } else { + for (size_t i = 0; i < rec.size(); ++i) { + rs_func_.push_back(std::make_shared()); + } + } - rv_func_.reserve(rec.size()); - if (FluidSystem::enableVaporizedOil()) { - const Opm::TableContainer& rvvdTables = tables.getRvvdTables(); - for (size_t i = 0; i < rec.size(); ++i) { - if (eqlmap.cells(i).empty()) - { - rv_func_.push_back(std::shared_ptr>()); - continue; - } - const int pvtIdx = regionPvtIdx_[i]; - if (!rec[i].wetGasInitConstantRv()) { - if (rvvdTables.size() <= 0) { - OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available."); - } + rv_func_.reserve(rec.size()); + if (FluidSystem::enableVaporizedOil()) { + const Opm::TableContainer& rvvdTables = tables.getRvvdTables(); + for (size_t i = 0; i < rec.size(); ++i) { + if (eqlmap.cells(i).empty()) + { + rv_func_.push_back(std::shared_ptr>()); + continue; + } + const int pvtIdx = regionPvtIdx_[i]; + if (!rec[i].wetGasInitConstantRv()) { + if (rvvdTables.size() <= 0) { + OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table not available."); + } - const Opm::RvvdTable& rvvdTable = rvvdTables.getTable(i); - std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); - std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); - rv_func_.push_back(std::make_shared>(pvtIdx, + const Opm::RvvdTable& rvvdTable = rvvdTables.getTable(i); + std::vector depthColumn = rvvdTable.getColumn("DEPTH").vectorCopy(); + std::vector rvColumn = rvvdTable.getColumn("RV").vectorCopy(); + rv_func_.push_back(std::make_shared>(pvtIdx, depthColumn , rvColumn)); - } else { - if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { - OPM_THROW(std::runtime_error, - "Cannot initialise: when no explicit RVVD table is given, \n" - "datum depth must be at the gas-oil-contact. " - "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); - } - const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); - const double T_contact = 273.15 + 20; // standard temperature for now - rv_func_.push_back(std::make_shared>(pvtIdx ,p_contact, T_contact)); - } - } - } else { - for (size_t i = 0; i < rec.size(); ++i) { - rv_func_.push_back(std::make_shared()); - } + } else { + if (rec[i].gasOilContactDepth() != rec[i].datumDepth()) { + OPM_THROW(std::runtime_error, + "Cannot initialise: when no explicit RVVD table is given, \n" + "datum depth must be at the gas-oil-contact. " + "In EQUIL region " << (i + 1) << " (counting from 1), this does not hold."); } - - // Compute pressures, saturations, rs and rv factors. - calcPressSatRsRv(eqlmap, rec, materialLawManager, G, grav); - - // Modify oil pressure in no-oil regions so that the pressures of present phases can - // be recovered from the oil pressure and capillary relations. + const double p_contact = rec[i].datumDepthPressure() + rec[i].gasOilContactCapillaryPressure(); + const double T_contact = 273.15 + 20; // standard temperature for now + rv_func_.push_back(std::make_shared>(pvtIdx ,p_contact, T_contact)); } + } + } else { + for (size_t i = 0; i < rec.size(); ++i) { + rv_func_.push_back(std::make_shared()); + } + } - typedef std::vector Vec; - typedef std::vector PVec; // One per phase. + // Compute pressures, saturations, rs and rv factors. + calcPressSatRsRv(eqlmap, rec, materialLawManager, grid, grav); - const PVec& press() const { return pp_; } - const PVec& saturation() const { return sat_; } - const Vec& rs() const { return rs_; } - const Vec& rv() const { return rv_; } + // Modify oil pressure in no-oil regions so that the pressures of present phases can + // be recovered from the oil pressure and capillary relations. + } - private: - typedef EquilReg EqReg; - std::vector< std::shared_ptr > rs_func_; - std::vector< std::shared_ptr > rv_func_; - std::vector regionPvtIdx_; - PVec pp_; - PVec sat_; - Vec rs_; - Vec rv_; - Vec swat_init_; + typedef std::vector Vec; + typedef std::vector PVec; // One per phase. - template - void setRegionPvtIdx(const Grid& G, const Opm::EclipseState& eclipseState, const RMap& reg) { + const PVec& press() const { return pp_; } + const PVec& saturation() const { return sat_; } + const Vec& rs() const { return rs_; } + const Vec& rv() const { return rv_; } - std::vector cellPvtRegionIdx; - extractPvtTableIndex(cellPvtRegionIdx, eclipseState, Opm::UgGridHelpers::numCells(G), Opm::UgGridHelpers::globalCell(G)); - for (const auto& r : reg.activeRegions()) { - const auto& cells = reg.cells(r); - const int cell = *(cells.begin()); - regionPvtIdx_[r] = cellPvtRegionIdx[cell]; - } - } +private: + typedef EquilReg EqReg; + std::vector< std::shared_ptr > rs_func_; + std::vector< std::shared_ptr > rv_func_; + std::vector regionPvtIdx_; + PVec pp_; + PVec sat_; + Vec rs_; + Vec rv_; + Vec swat_init_; - template - void - calcPressSatRsRv(const RMap& reg , - const std::vector< Opm::EquilRecord >& rec , - MaterialLawManager& materialLawManager, - const Grid& G , - const double grav) - { - for (const auto& r : reg.activeRegions()) { - const auto& cells = reg.cells(r); - if (cells.empty()) - { - Opm::OpmLog::warning("Equilibration region " + std::to_string(r + 1) - + " has no active cells"); - continue; - } + template + void setRegionPvtIdx(const Grid& grid, const Opm::EclipseState& eclipseState, const RMap& reg) { - const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]); - - PVec pressures = phasePressures(G, eqreg, cells, grav); - const std::vector& temp = temperature(G, eqreg, cells); - const PVec sat = phaseSaturations(G, eqreg, cells, materialLawManager, swat_init_, pressures); + std::vector cellPvtRegionIdx; + extractPvtTableIndex(cellPvtRegionIdx, eclipseState, Opm::UgGridHelpers::numCells(grid), Opm::UgGridHelpers::globalCell(grid)); + for (const auto& r : reg.activeRegions()) { + const auto& cells = reg.cells(r); + const int cell = *(cells.begin()); + regionPvtIdx_[r] = cellPvtRegionIdx[cell]; + } + } - const int np = FluidSystem::numPhases; - for (int p = 0; p < np; ++p) { - copyFromRegion(pressures[p], cells, pp_[p]); - copyFromRegion(sat[p], cells, sat_[p]); - } - const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - if (oil && gas) { - const int oilpos = FluidSystem::oilPhaseIdx; - const int gaspos = FluidSystem::gasPhaseIdx; - const Vec rs_vals = computeRs(G, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]); - const Vec rv_vals = computeRs(G, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]); - copyFromRegion(rs_vals, cells, rs_); - copyFromRegion(rv_vals, cells, rv_); - } - } - } + template + void + calcPressSatRsRv(const RMap& reg , + const std::vector< Opm::EquilRecord >& rec , + MaterialLawManager& materialLawManager, + const Grid& grid , + const double grav) + { + for (const auto& r : reg.activeRegions()) { + const auto& cells = reg.cells(r); + if (cells.empty()) + { + Opm::OpmLog::warning("Equilibration region " + std::to_string(r + 1) + + " has no active cells"); + continue; + } - template - void copyFromRegion(const Vec& source, - const CellRangeType& cells, - Vec& destination) - { - auto s = source.begin(); - auto c = cells.begin(); - const auto e = cells.end(); - for (; c != e; ++c, ++s) { - destination[*c] = *s; - } - } + const EqReg eqreg(rec[r], rs_func_[r], rv_func_[r], regionPvtIdx_[r]); - }; - } // namespace DeckDependent - } // namespace EQUIL + PVec pressures = phasePressures(grid, eqreg, cells, grav); + const std::vector& temp = temperature(grid, eqreg, cells); + const PVec sat = phaseSaturations(grid, eqreg, cells, materialLawManager, swat_init_, pressures); + + const int np = FluidSystem::numPhases; + for (int p = 0; p < np; ++p) { + copyFromRegion(pressures[p], cells, pp_[p]); + copyFromRegion(sat[p], cells, sat_[p]); + } + const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); + const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); + if (oil && gas) { + const int oilpos = FluidSystem::oilPhaseIdx; + const int gaspos = FluidSystem::gasPhaseIdx; + const Vec rs_vals = computeRs(grid, cells, pressures[oilpos], temp, *(rs_func_[r]), sat[gaspos]); + const Vec rv_vals = computeRs(grid, cells, pressures[gaspos], temp, *(rv_func_[r]), sat[oilpos]); + copyFromRegion(rs_vals, cells, rs_); + copyFromRegion(rv_vals, cells, rv_); + } + } + } + + template + void copyFromRegion(const Vec& source, + const CellRangeType& cells, + Vec& destination) + { + auto s = source.begin(); + auto c = cells.begin(); + const auto e = cells.end(); + for (; c != e; ++c, ++s) { + destination[*c] = *s; + } + } + +}; +} // namespace DeckDependent +} // namespace EQUIL } // namespace Opm -#include "initStateEquil_impl.hpp" - #endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED diff --git a/ebos/equil/initStateEquil_impl.hpp b/ebos/equil/initStateEquil_impl.hpp deleted file mode 100644 index 8e24021b2..000000000 --- a/ebos/equil/initStateEquil_impl.hpp +++ /dev/null @@ -1,796 +0,0 @@ -// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- -// vi: set et ts=4 sw=4 sts=4: -/* - This file is part of the Open Porous Media project (OPM). - - OPM is free software: you can redistribute it and/or modify - it under the terms of the GNU General Public License as published by - the Free Software Foundation, either version 3 of the License, or - (at your option) any later version. - - OPM is distributed in the hope that it will be useful, - but WITHOUT ANY WARRANTY; without even the implied warranty of - MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - GNU General Public License for more details. - - You should have received a copy of the GNU General Public License - along with OPM. If not, see . - - Consult the COPYING file in the top-level source directory of this - module for the precise wording of the license and the list of - copyright holders. -*/ -/** - * \file - * - * \brief Routines that actually solve the ODEs that emerge from the hydrostatic - * equilibrium problem - */ -#ifndef EWOMS_INITSTATEEQUIL_IMPL_HEADER_INCLUDED -#define EWOMS_INITSTATEEQUIL_IMPL_HEADER_INCLUDED - -#include - -#include - -#include -#include -#include -#include - -namespace Ewoms -{ - namespace Details { - template - class RK4IVP { - public: - RK4IVP(const RHS& f , - const std::array& span, - const double y0 , - const int N ) - : N_(N) - , span_(span) - { - const double h = stepsize(); - const double h2 = h / 2; - const double h6 = h / 6; - - y_.reserve(N + 1); - f_.reserve(N + 1); - - y_.push_back(y0); - f_.push_back(f(span_[0], y0)); - - for (int i = 0; i < N; ++i) { - const double x = span_[0] + i*h; - const double y = y_.back(); - - const double k1 = f_[i]; - const double k2 = f(x + h2, y + h2*k1); - const double k3 = f(x + h2, y + h2*k2); - const double k4 = f(x + h , y + h*k3); - - y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4)); - f_.push_back(f(x + h, y_.back())); - } - - assert (y_.size() == std::vector::size_type(N + 1)); - } - - double - operator()(const double x) const - { - // Dense output (O(h**3)) according to Shampine - // (Hermite interpolation) - const double h = stepsize(); - int i = (x - span_[0]) / h; - const double t = (x - (span_[0] + i*h)) / h; - - // Crude handling of evaluation point outside "span_"; - if (i < 0) { i = 0; } - if (N_ <= i) { i = N_ - 1; } - - const double y0 = y_[i], y1 = y_[i + 1]; - const double f0 = f_[i], f1 = f_[i + 1]; - - double u = (1 - 2*t) * (y1 - y0); - u += h * ((t - 1)*f0 + t*f1); - u *= t * (t - 1); - u += (1 - t)*y0 + t*y1; - - return u; - } - - private: - int N_; - std::array span_; - std::vector y_; - std::vector f_; - - double - stepsize() const { return (span_[1] - span_[0]) / N_; } - }; - - namespace PhasePressODE { - template - class Water { - public: - Water(const double temp, - const int pvtRegionIdx, - const double norm_grav) - : temp_(temp) - , pvtRegionIdx_(pvtRegionIdx) - , g_(norm_grav) - { - } - - double - operator()(const double /* depth */, - const double press) const - { - return this->density(press) * g_; - } - - private: - const double temp_; - const int pvtRegionIdx_; - const double g_; - - double - density(const double press) const - { - double rho = FluidSystem::waterPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - rho *= FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, pvtRegionIdx_); - return rho; - } - }; - - template - class Oil { - public: - Oil(const double temp, - const RS& rs, - const int pvtRegionIdx, - const double norm_grav) - : temp_(temp) - , rs_(rs) - , pvtRegionIdx_(pvtRegionIdx) - , g_(norm_grav) - { - } - - double - operator()(const double depth, - const double press) const - { - return this->density(depth, press) * g_; - } - - private: - const double temp_; - const RS& rs_; - const int pvtRegionIdx_; - const double g_; - - double - density(const double depth, - const double press) const - { - double rs = rs_(depth, press, temp_); - double bOil = 0.0; - if ( !FluidSystem::enableDissolvedGas() || rs >= FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp_, press) ) { - bOil = FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { - bOil = FluidSystem::oilPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rs); - } - double rho = bOil * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - if (FluidSystem::enableDissolvedGas()) { - rho += rs * bOil * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - } - - return rho; - } - }; - - template - class Gas { - public: - Gas(const double temp, - const RV& rv, - const int pvtRegionIdx, - const double norm_grav) - : temp_(temp) - , rv_(rv) - , pvtRegionIdx_(pvtRegionIdx) - , g_(norm_grav) - { - } - - double - operator()(const double depth, - const double press) const - { - return this->density(depth, press) * g_; - } - - private: - const double temp_; - const RV& rv_; - const int pvtRegionIdx_; - const double g_; - - double - density(const double depth, - const double press) const - { - double rv = rv_(depth, press, temp_); - double bGas = 0.0; - if ( !FluidSystem::enableVaporizedOil() || rv >= FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp_, press) ) { - bGas = FluidSystem::gasPvt().saturatedInverseFormationVolumeFactor(pvtRegionIdx_, temp_, press); - } else { - bGas = FluidSystem::gasPvt().inverseFormationVolumeFactor(pvtRegionIdx_, temp_, press, rv); - } - double rho = bGas * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, pvtRegionIdx_); - if (FluidSystem::enableVaporizedOil()) { - rho += rv * bGas * FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, pvtRegionIdx_); - } - - return rho; - } - }; - } // namespace PhasePressODE - - - namespace PhasePressure { - template - void - assign(const Grid& G , - const std::array& f , - const double split, - const CellRange& cells, - std::vector& p ) - { - - enum { up = 0, down = 1 }; - - std::vector::size_type c = 0; - for (typename CellRange::const_iterator - ci = cells.begin(), ce = cells.end(); - ci != ce; ++ci, ++c) - { - assert (c < p.size()); - - const double z = Opm::UgGridHelpers::cellCenterDepth(G, *ci); - p[c] = (z < split) ? f[up](z) : f[down](z); - } - } - - template < class FluidSystem, - class Grid, - class Region, - class CellRange> - void - water(const Grid& G , - const Region& reg , - const std::array& span , - const double grav , - double& po_woc, - const CellRange& cells , - std::vector& press ) - { - using PhasePressODE::Water; - typedef Water ODE; - - const double T = 273.15 + 20; // standard temperature for now - ODE drho(T, reg.pvtIdx() , grav); - - double z0; - double p0; - if (reg.datum() > reg.zwoc()) {//Datum in water zone - z0 = reg.datum(); - p0 = reg.pressure(); - } else { - z0 = reg.zwoc(); - p0 = po_woc - reg.pcow_woc(); // Water pressure at contact - } - - std::array up = {{ z0, span[0] }}; - std::array down = {{ z0, span[1] }}; - - typedef Details::RK4IVP WPress; - std::array wpress = { - { - WPress(drho, up , p0, 2000) - , - WPress(drho, down, p0, 2000) - } - }; - - assign(G, wpress, z0, cells, press); - - if (reg.datum() > reg.zwoc()) { - // Return oil pressure at contact - po_woc = wpress[0](reg.zwoc()) + reg.pcow_woc(); - } - } - - template - void - oil(const Grid& G , - const Region& reg , - const std::array& span , - const double grav , - const CellRange& cells , - std::vector& press , - double& po_woc, - double& po_goc) - { - using PhasePressODE::Oil; - typedef Oil ODE; - - const double T = 273.15 + 20; // standard temperature for now - ODE drho(T, reg.dissolutionCalculator(), - reg.pvtIdx(), grav); - - double z0; - double p0; - if (reg.datum() > reg.zwoc()) {//Datum in water zone, po_woc given - z0 = reg.zwoc(); - p0 = po_woc; - } else if (reg.datum() < reg.zgoc()) {//Datum in gas zone, po_goc given - z0 = reg.zgoc(); - p0 = po_goc; - } else { //Datum in oil zone - z0 = reg.datum(); - p0 = reg.pressure(); - } - - std::array up = {{ z0, span[0] }}; - std::array down = {{ z0, span[1] }}; - - typedef Details::RK4IVP OPress; - std::array opress = { - { - OPress(drho, up , p0, 2000) - , - OPress(drho, down, p0, 2000) - } - }; - - assign(G, opress, z0, cells, press); - - const double woc = reg.zwoc(); - if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum - else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum - else { po_woc = p0; } // WOC *at* datum - - const double goc = reg.zgoc(); - if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum - else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum - else { po_goc = p0; } // GOC *at* datum - } - - template - void - gas(const Grid& G , - const Region& reg , - const std::array& span , - const double grav , - double& po_goc, - const CellRange& cells , - std::vector& press ) - { - using PhasePressODE::Gas; - typedef Gas ODE; - - const double T = 273.15 + 20; // standard temperature for now - ODE drho(T, reg.evaporationCalculator(), - reg.pvtIdx(), grav); - - double z0; - double p0; - if (reg.datum() < reg.zgoc()) {//Datum in gas zone - z0 = reg.datum(); - p0 = reg.pressure(); - } else { - z0 = reg.zgoc(); - p0 = po_goc + reg.pcgo_goc(); // Gas pressure at contact - } - - std::array up = {{ z0, span[0] }}; - std::array down = {{ z0, span[1] }}; - - typedef Details::RK4IVP GPress; - std::array gpress = { - { - GPress(drho, up , p0, 2000) - , - GPress(drho, down, p0, 2000) - } - }; - - assign(G, gpress, z0, cells, press); - - if (reg.datum() < reg.zgoc()) { - // Return oil pressure at contact - po_goc = gpress[1](reg.zgoc()) - reg.pcgo_goc(); - } - } - } // namespace PhasePressure - - template - void - equilibrateOWG(const Grid& G, - const Region& reg, - const double grav, - const std::array& span, - const CellRange& cells, - std::vector< std::vector >& press) - { - const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); - const bool oil = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx); - const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - const int oilpos = FluidSystem::oilPhaseIdx; - const int waterpos = FluidSystem::waterPhaseIdx; - const int gaspos = FluidSystem::gasPhaseIdx; - - if (reg.datum() > reg.zwoc()) { // Datum in water zone - double po_woc = -1; - double po_goc = -1; - - if (water) { - PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ waterpos ]); - } - - if (oil) { - PhasePressure::oil(G, reg, span, grav, cells, - press[ oilpos ], po_woc, po_goc); - } - - if (gas) { - PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gaspos ]); - } - } else if (reg.datum() < reg.zgoc()) { // Datum in gas zone - double po_woc = -1; - double po_goc = -1; - - if (gas) { - PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gaspos ]); - } - - if (oil) { - PhasePressure::oil(G, reg, span, grav, cells, - press[ oilpos ], po_woc, po_goc); - } - - if (water) { - PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ waterpos ]); - } - } else { // Datum in oil zone - double po_woc = -1; - double po_goc = -1; - - if (oil) { - PhasePressure::oil(G, reg, span, grav, cells, - press[ oilpos ], po_woc, po_goc); - } - - if (water) { - PhasePressure::water(G, reg, span, grav, po_woc, - cells, press[ waterpos ]); - } - - if (gas) { - PhasePressure::gas(G, reg, span, grav, po_goc, - cells, press[ gaspos ]); - } - } - } - } // namespace Details - - - namespace EQUIL { - - - template - std::vector< std::vector > - phasePressures(const Grid& G, - const Region& reg, - const CellRange& cells, - const double grav) - { - std::array span = - {{ std::numeric_limits::max() , - -std::numeric_limits::max() }}; // Symm. about 0. - - int ncell = 0; - { - // This code is only supported in three space dimensions - assert (Opm::UgGridHelpers::dimensions(G) == 3); - - const int nd = Opm::UgGridHelpers::dimensions(G); - - // Define vertical span as - // - // [minimum(node depth(cells)), maximum(node depth(cells))] - // - // Note: We use a sledgehammer approach--looping all - // the nodes of all the faces of all the 'cells'--to - // compute those bounds. This necessarily entails - // visiting some nodes (and faces) multiple times. - // - // Note: The implementation of 'RK4IVP<>' implicitly - // imposes the requirement that cell centroids are all - // within this vertical span. That requirement is not - // checked. - auto cell2Faces = Opm::UgGridHelpers::cell2Faces(G); - auto faceVertices = Opm::UgGridHelpers::face2Vertices(G); - - for (typename CellRange::const_iterator - ci = cells.begin(), ce = cells.end(); - ci != ce; ++ci, ++ncell) - { - for (auto fi=cell2Faces[*ci].begin(), - fe=cell2Faces[*ci].end(); - fi != fe; - ++fi) - { - for (auto i = faceVertices[*fi].begin(), e = faceVertices[*fi].end(); - i != e; ++i) - { - const double z = Opm::UgGridHelpers::vertexCoordinates(G, *i)[nd-1]; - - if (z < span[0]) { span[0] = z; } - if (z > span[1]) { span[1] = z; } - } - } - } - } - const int np = FluidSystem::numPhases; //reg.phaseUsage().num_phases; - - typedef std::vector pval; - std::vector press(np, pval(ncell, 0.0)); - - const double zwoc = reg.zwoc (); - const double zgoc = reg.zgoc (); - - // make sure goc and woc is within the span for the phase pressure calculation - span[0] = std::min(span[0],zgoc); - span[1] = std::max(span[1],zwoc); - - Details::equilibrateOWG(G, reg, grav, span, cells, press); - - return press; - } - - template - std::vector - temperature(const Grid& /* G */, - const Region& /* reg */, - const CellRange& cells) - { - // use the standard temperature for everything for now - return std::vector(cells.size(), 273.15 + 20.0); - } - - template - std::vector< std::vector > - phaseSaturations(const Grid& G, - const Region& reg, - const CellRange& cells, - MaterialLawManager& materialLawManager, - const std::vector swat_init, - std::vector< std::vector >& phase_pressures) - { - if (!FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) { - OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases."); - } - - std::vector< std::vector > phase_saturations = phase_pressures; // Just to get the right size. - - // Adjust oil pressure according to gas saturation and cap pressure - typedef Opm::SimpleModularFluidState SatOnlyFluidState; - - SatOnlyFluidState fluidState; - typedef typename MaterialLawManager::MaterialLaw MaterialLaw; - - const bool water = FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx); - const bool gas = FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx); - const int oilpos = FluidSystem::oilPhaseIdx; - const int waterpos = FluidSystem::waterPhaseIdx; - const int gaspos = FluidSystem::gasPhaseIdx; - std::vector::size_type local_index = 0; - for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) { - const int cell = *ci; - const auto& scaledDrainageInfo = - materialLawManager.oilWaterScaledEpsInfoDrainage(cell); - const auto& matParams = materialLawManager.materialLawParams(cell); - - // Find saturations from pressure differences by - // inverting capillary pressure functions. - double sw = 0.0; - if (water) { - if (isConstPc(materialLawManager,FluidSystem::waterPhaseIdx, cell)){ - const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(G, - cell); - sw = satFromDepth(materialLawManager,cellDepth,reg.zwoc(),waterpos,cell,false); - phase_saturations[waterpos][local_index] = sw; - } - else{ - const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index]; - if (swat_init.empty()) { // Invert Pc to find sw - sw = satFromPc(materialLawManager, waterpos, cell, pcov); - phase_saturations[waterpos][local_index] = sw; - } else { // Scale Pc to reflect imposed sw - sw = swat_init[cell]; - sw = materialLawManager.applySwatinit(cell, pcov, sw); - phase_saturations[waterpos][local_index] = sw; - } - } - } - double sg = 0.0; - if (gas) { - if (isConstPc(materialLawManager,FluidSystem::gasPhaseIdx,cell)){ - const double cellDepth = Opm::UgGridHelpers::cellCenterDepth(G, - cell); - sg = satFromDepth(materialLawManager,cellDepth,reg.zgoc(),gaspos,cell,true); - phase_saturations[gaspos][local_index] = sg; - } - else{ - // Note that pcog is defined to be (pg - po), not (po - pg). - const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index]; - const double increasing = true; // pcog(sg) expected to be increasing function - sg = satFromPc(materialLawManager, gaspos, cell, pcog, increasing); - phase_saturations[gaspos][local_index] = sg; - } - } - if (gas && water && (sg + sw > 1.0)) { - // Overlapping gas-oil and oil-water transition - // zones can lead to unphysical saturations when - // treated as above. Must recalculate using gas-water - // capillary pressure. - const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index]; - if (! swat_init.empty()) { - // Re-scale Pc to reflect imposed sw for vanishing oil phase. - // This seems consistent with ecl, and fails to honour - // swat_init in case of non-trivial gas-oil cap pressure. - sw = materialLawManager.applySwatinit(cell, pcgw, sw); - } - sw = satFromSumOfPcs(materialLawManager, waterpos, gaspos, cell, pcgw); - sg = 1.0 - sw; - phase_saturations[waterpos][local_index] = sw; - phase_saturations[gaspos][local_index] = sg; - if ( water ) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, sw); - } - else { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); - } - fluidState.setSaturation(FluidSystem::oilPhaseIdx, 1.0 - sw - sg); - fluidState.setSaturation(FluidSystem::gasPhaseIdx, sg); - - double pC[/*numPhases=*/3] = { 0.0, 0.0, 0.0 }; - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; - phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas; - } - phase_saturations[oilpos][local_index] = 1.0 - sw - sg; - - // Adjust phase pressures for max and min saturation ... - double threshold_sat = 1.0e-6; - - double so = 1.0; - double pC[FluidSystem::numPhases] = { 0.0, 0.0, 0.0 }; - if (water) { - double swu = scaledDrainageInfo.Swu; - fluidState.setSaturation(FluidSystem::waterPhaseIdx, swu); - so -= swu; - } - if (gas) { - double sgu = scaledDrainageInfo.Sgu; - fluidState.setSaturation(FluidSystem::gasPhaseIdx, sgu); - so-= sgu; - } - fluidState.setSaturation(FluidSystem::oilPhaseIdx, so); - - if (water && sw > scaledDrainageInfo.Swu-threshold_sat ) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swu); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; - phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pcWat; - } else if (gas && sg > scaledDrainageInfo.Sgu-threshold_sat) { - fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgu); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; - phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pcGas; - } - if (gas && sg < scaledDrainageInfo.Sgl+threshold_sat) { - fluidState.setSaturation(FluidSystem::gasPhaseIdx, scaledDrainageInfo.Sgl); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcGas = pC[FluidSystem::oilPhaseIdx] + pC[FluidSystem::gasPhaseIdx]; - phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pcGas; - } - if (water && sw < scaledDrainageInfo.Swl+threshold_sat) { - fluidState.setSaturation(FluidSystem::waterPhaseIdx, scaledDrainageInfo.Swl); - MaterialLaw::capillaryPressures(pC, matParams, fluidState); - double pcWat = pC[FluidSystem::oilPhaseIdx] - pC[FluidSystem::waterPhaseIdx]; - phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pcWat; - } - } - return phase_saturations; - } - - - /** - * Compute initial Rs values. - * - * \tparam CellRangeType Type of cell range that demarcates the - * cells pertaining to the current - * equilibration region. Must implement - * methods begin() and end() to bound the range - * as well as provide an inner type, - * const_iterator, to traverse the range. - * - * \param[in] grid Grid. - * \param[in] cells Range that spans the cells of the current - * equilibration region. - * \param[in] oil_pressure Oil pressure for each cell in range. - * \param[in] temperature Temperature for each cell in range. - * \param[in] rs_func Rs as function of pressure and depth. - * \return Rs values, one for each cell in the 'cells' range. - */ - template - std::vector computeRs(const Grid& grid, - const CellRangeType& cells, - const std::vector oil_pressure, - const std::vector& temperature, - const Miscibility::RsFunction& rs_func, - const std::vector gas_saturation) - { - assert(Opm::UgGridHelpers::dimensions(grid) == 3); - std::vector rs(cells.size()); - int count = 0; - for (auto it = cells.begin(); it != cells.end(); ++it, ++count) { - const double depth = Opm::UgGridHelpers::cellCenterDepth(grid, *it); - rs[count] = rs_func(depth, oil_pressure[count], temperature[count], gas_saturation[count]); - } - return rs; - } - - } // namespace Equil - - -} // namespace Opm - -#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED