diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 3aaa601c..fb37a7f8 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -233,6 +233,12 @@ namespace Opm const double pcow, double & swat); + // return a reference to the "raw" PVT fluid object for a phase. + const PvtInterface& pvt(int phaseIdx) const + { + return pvt_.pvt(phaseIdx); + } + private: int getTableIndex_(const int* pvtTableIdx, int cellIdx) const { diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp index 45caa12f..c1778094 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.hpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.hpp @@ -117,6 +117,12 @@ namespace Opm double* output_R, double* output_dRdp) const; + // return a reference to the "raw" PVT fluid object for a phase. + const PvtInterface& pvt(int phaseIdx) const + { + return *props_[phaseIdx]; + } + private: // Disabling copying (just to avoid surprises, since we use shared_ptr). BlackoilPvtProperties(const BlackoilPvtProperties&); diff --git a/opm/core/utility/thresholdPressures.hpp b/opm/core/utility/thresholdPressures.hpp index 8fda8cdc..1ae5e899 100644 --- a/opm/core/utility/thresholdPressures.hpp +++ b/opm/core/utility/thresholdPressures.hpp @@ -20,10 +20,14 @@ #ifndef OPM_THRESHOLDPRESSURES_HEADER_INCLUDED #define OPM_THRESHOLDPRESSURES_HEADER_INCLUDED +#include +#include + #include #include #include -#include +#include +#include namespace Opm @@ -48,8 +52,15 @@ namespace Opm template - std::vector thresholdPressures(const ParseMode& parseMode ,EclipseStateConstPtr eclipseState, const Grid& grid) + std::vector thresholdPressures(const DeckConstPtr& deck, + EclipseStateConstPtr eclipseState, + const Grid& grid, + const BlackoilState& initialState, + const BlackoilPropertiesFromDeck& props, + const double gravity) { + + const PhaseUsage& pu = props.phaseUsage(); SimulationConfigConstPtr simulationConfig = eclipseState->getSimulationConfig(); std::vector thpres_vals; if (simulationConfig->hasThresholdPressure()) { @@ -57,7 +68,172 @@ namespace Opm std::shared_ptr> eqlnum = eclipseState->getIntGridProperty("EQLNUM"); auto eqlnumData = eqlnum->getData(); - // Set values for each face. + int numPhases = initialState.numPhases(); + int numCells = UgGridHelpers::numCells(grid); + int numPvtRegions = deck->getKeyword("TABDIMS")->getRecord(0)->getItem("NTPVT")->getInt(0); + + // retrieve the surface densities + std::vector > surfaceDensity(numPvtRegions); + Opm::DeckKeywordConstPtr densityKw = deck->getKeyword("DENSITY"); + for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) { + surfaceDensity[regionIdx].resize(numPhases); + + if (pu.phase_used[BlackoilPhases::Aqua]) { + int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + surfaceDensity[regionIdx][wpos] = + densityKw->getRecord(regionIdx)->getItem("WATER")->getSIDouble(0); + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + surfaceDensity[regionIdx][opos] = + densityKw->getRecord(regionIdx)->getItem("OIL")->getSIDouble(0); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + surfaceDensity[regionIdx][gpos] = + densityKw->getRecord(regionIdx)->getItem("GAS")->getSIDouble(0); + } + } + + // retrieve the PVT region of each cell. note that we need c++ instead of + // Fortran indices. + std::vector pvtRegion = eclipseState->getIntGridProperty("PVTNUM")->getData(); + for (auto it = pvtRegion.begin(); it != pvtRegion.end(); ++ it) { + *it = std::max(0, *it - 1); + } + + // compute the initial "phase presence" of each cell (required to calculate + // the inverse formation volume factors + std::vector cond(numCells); + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + if (pu.phase_used[BlackoilPhases::Aqua] && initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]] > 0.0) { + cond[cellIdx].setFreeWater(); + } + + if (pu.phase_used[BlackoilPhases::Liquid] && initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]] > 0.0) { + cond[cellIdx].setFreeOil(); + } + + if (pu.phase_used[BlackoilPhases::Vapour] && initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]] > 0.0) { + cond[cellIdx].setFreeGas(); + } + } + + // calculate the initial fluid densities for the gravity correction. + std::vector b(numCells); + + std::vector> rho(numPhases); + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + rho[phaseIdx].resize(numCells); + } + + // compute the capillary pressures of the active phases + std::vector capPress(numCells*numPhases); + std::vector cellIdxArray(numCells); + for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { + cellIdxArray[cellIdx] = cellIdx; + } + props.capPress(numCells, initialState.saturation().data(), cellIdxArray.data(), capPress.data(), NULL); + + // compute the absolute pressure of each active phase: for some reason, E100 + // defines the capillary pressure for the water phase as p_o - p_w while it + // uses p_g - p_o for the gas phase. (it would be more consistent to use the + // oil pressure as reference for both the other phases.) probably this is + // done to always have a positive number for the capillary pressure (as long + // as the medium is hydrophilic) + std::vector > phasePressure(numPhases); + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + phasePressure[phaseIdx].resize(numCells); + } + + for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { + assert(pu.phase_used[BlackoilPhases::Liquid]); // we currently hard-code the oil phase as the reference phase! + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx]; + + if (pu.phase_used[BlackoilPhases::Aqua]) { + int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + phasePressure[wpos][cellIdx] = initialState.pressure()[cellIdx] + (capPress[cellIdx*numPhases + opos] - capPress[cellIdx*numPhases + wpos]); + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + phasePressure[gpos][cellIdx] = initialState.pressure()[cellIdx] + (capPress[cellIdx*numPhases + gpos] - capPress[cellIdx*numPhases + opos]); + } + } + + // calculate the inverse formation volume factors for the active phases and each cell + if (pu.phase_used[BlackoilPhases::Aqua]) { + std::vector dummy(numCells*BlackoilPhases::MaxNumPhases); + int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const PvtInterface& pvtw = props.pvt(wpos); + pvtw.b(numCells, + pvtRegion.data(), + phasePressure[wpos].data(), + initialState.temperature().data(), + initialState.gasoilratio().data(), + cond.data(), + b.data(), + dummy.data(), + dummy.data()); + for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { + rho[wpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][wpos]*b[cellIdx]; + } + } + + if (pu.phase_used[BlackoilPhases::Liquid]) { + std::vector dummy(numCells*BlackoilPhases::MaxNumPhases); + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const PvtInterface& pvto = props.pvt(opos); + pvto.b(numCells, + pvtRegion.data(), + phasePressure[opos].data(), + initialState.temperature().data(), + initialState.gasoilratio().data(), + cond.data(), + b.data(), + dummy.data(), + dummy.data()); + for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { + rho[opos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][opos]*b[cellIdx]; + + if (pu.phase_used[BlackoilPhases::Vapour]) { + int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + rho[opos][cellIdx] += + surfaceDensity[pvtRegion[cellIdx]][gpos]*initialState.gasoilratio()[cellIdx]*b[cellIdx]; + } + } + } + + if (pu.phase_used[BlackoilPhases::Vapour]) { + std::vector dummy(numCells*BlackoilPhases::MaxNumPhases); + int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const PvtInterface& pvtg = props.pvt(gpos); + pvtg.b(numCells, + pvtRegion.data(), + phasePressure[gpos].data(), + initialState.temperature().data(), + initialState.gasoilratio().data(), + cond.data(), + b.data(), + dummy.data(), + dummy.data()); + for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) { + rho[gpos][cellIdx] = surfaceDensity[pvtRegion[cellIdx]][gpos]*b[cellIdx]; + + if (pu.phase_used[BlackoilPhases::Liquid]) { + int opos = pu.phase_pos[BlackoilPhases::Liquid]; + rho[gpos][cellIdx] += + surfaceDensity[pvtRegion[cellIdx]][opos]*initialState.rv()[cellIdx]*b[cellIdx]; + } + } + } + + // Calculate the maximum pressure potential difference between all PVT region + // transitions of the initial solution. + std::map, double> maxDp; const int num_faces = UgGridHelpers::numFaces(grid); thpres_vals.resize(num_faces, 0.0); const int* gc = UgGridHelpers::globalCell(grid); @@ -74,12 +250,64 @@ namespace Opm const int eq1 = eqlnumData[gc1]; const int eq2 = eqlnumData[gc2]; + if (eq1 == eq2) { + // not an equilibration region boundary. skip this. + continue; + } + + // update the maximum pressure potential difference between the two + // regions + const auto barrierId = std::make_pair(eq1, eq2); + if (maxDp.count(barrierId) == 0) { + maxDp[barrierId] = 0.0; + } + + for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + double z1 = UgGridHelpers::cellCenterDepth(grid, c1); + double z2 = UgGridHelpers::cellCenterDepth(grid, c2); + double zAvg = (z1 + z2)/2; // average depth + + double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2; + + double s1 = initialState.saturation()[numPhases*c1 + phaseIdx]; + double s2 = initialState.saturation()[numPhases*c2 + phaseIdx]; + + // compute gravity corrected pressure potentials at the average depth + double p1 = phasePressure[phaseIdx][c1]; + double p2 = phasePressure[phaseIdx][c2]; + + p1 += rhoAvg*gravity*(z1 - zAvg); + p2 += rhoAvg*gravity*(z2 - zAvg); + + if ((p1 > p2 && s1 > 0.0) || (p2 > p1 && s2 > 0.0)) + maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2)); + } + } + + // Set threshold pressure values for each cell face. + thpres_vals.resize(num_faces, 0.0); + for (int face = 0; face < num_faces; ++face) { + const int c1 = fc(face, 0); + const int c2 = fc(face, 1); + if (c1 < 0 || c2 < 0) { + // Boundary face, skip it. + continue; + } + const int gc1 = (gc == 0) ? c1 : gc[c1]; + const int gc2 = (gc == 0) ? c2 : gc[c2]; + const int eq1 = eqlnumData[gc1]; + const int eq2 = eqlnumData[gc2]; + if (thresholdPressure->hasRegionBarrier(eq1,eq2)) { if (thresholdPressure->hasThresholdPressure(eq1,eq2)) { thpres_vals[face] = thresholdPressure->getThresholdPressure(eq1,eq2); - } else { - std::string msg = "Initializing the THPRES pressure values from the initial state is not supported - using 0.0"; - parseMode.handleError( ParseMode::UNSUPPORTED_INITIAL_THPRES , msg ); + } + else { + // set the threshold pressure for faces of PVT regions where the third item + // has been defaulted to the maximum pressure potential difference between + // these regions + const auto barrierId = std::make_pair(eq1, eq2); + thpres_vals[face] = maxDp[barrierId]; } } }