From 547b701a574094c42eac6ef91a348cedff1fc1ba Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Tue, 20 Oct 2015 13:37:44 +0200 Subject: [PATCH] implement determining the threshold pressure from the initial condition This needs to be done if a equilibration region transition is mentioned by the THPRES keyword, but no value is given for this record in the third item. (it seems that this is used quite frequently.) Also, the approach taken by this patch also does not collide with the restart machinery as far as I can see. This is because the initial condition is applied by the simulator before the state at the restart time is loaded. (I interpreted the code that way, but I could be wrong, could anyone verify this?) since it is pretty elaborate to calculate initial condition, this patch is pretty messy. I also do not know if Eclipse does include capillary pressure in this calculation or not (this patch does). Huge kudos go to [at]totto82 for reviewing, testing and debugging this. --- opm/core/props/BlackoilPropertiesFromDeck.hpp | 6 + opm/core/props/pvt/BlackoilPvtProperties.hpp | 6 + opm/core/utility/thresholdPressures.hpp | 240 +++++++++++++++++- 3 files changed, 246 insertions(+), 6 deletions(-) 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]; } } }