20 #ifndef OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
21 #define OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
23 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
24 #include <opm/core/props/BlackoilPhases.hpp>
27 #include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
28 #include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
29 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
30 #include <opm/parser/eclipse/Deck/Deck.hpp>
31 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
49 const EclipseState& eclipseState,
51 const BlackoilState& initialState,
52 const BlackoilPropertiesFromDeck& props,
56 const PhaseUsage& pu = props.phaseUsage();
58 const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty(
"EQLNUM");
59 const auto& eqlnumData = eqlnum.getData();
61 const int numPhases = initialState.numPhases();
62 const int numCells = UgGridHelpers::numCells(grid);
63 const int numPvtRegions = deck.getKeyword(
"TABDIMS").getRecord(0).getItem(
"NTPVT").get<
int >(0);
66 std::vector<double> minSat(numPhases*numCells);
67 std::vector<double> maxSat(numPhases*numCells);
68 std::vector<int> allCells(numCells);
69 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
70 allCells[cellIdx] = cellIdx;
72 props.satRange(numCells, allCells.data(), minSat.data(), maxSat.data());
75 std::vector<std::vector<double> > surfaceDensity(numPvtRegions);
76 const auto& densityKw = deck.getKeyword(
"DENSITY");
77 for (
int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) {
78 surfaceDensity[regionIdx].resize(numPhases);
80 if (pu.phase_used[BlackoilPhases::Aqua]) {
81 const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
82 surfaceDensity[regionIdx][wpos] =
83 densityKw.getRecord(regionIdx).getItem(
"WATER").getSIDouble(0);
86 if (pu.phase_used[BlackoilPhases::Liquid]) {
87 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
88 surfaceDensity[regionIdx][opos] =
89 densityKw.getRecord(regionIdx).getItem(
"OIL").getSIDouble(0);
92 if (pu.phase_used[BlackoilPhases::Vapour]) {
93 const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
94 surfaceDensity[regionIdx][gpos] =
95 densityKw.getRecord(regionIdx).getItem(
"GAS").getSIDouble(0);
101 const int* gc = UgGridHelpers::globalCell(grid);
102 std::vector<int> pvtRegion(numCells);
103 const auto& cartPvtRegion = eclipseState.get3DProperties().getIntGridProperty(
"PVTNUM").getData();
104 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
105 const int cartCellIdx = gc ? gc[cellIdx] : cellIdx;
106 pvtRegion[cellIdx] = std::max(0, cartPvtRegion[cartCellIdx] - 1);
111 std::vector<PhasePresence> cond(numCells);
112 for (
int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
113 if (pu.phase_used[BlackoilPhases::Aqua]) {
114 const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
116 cond[cellIdx].setFreeWater();
120 if (pu.phase_used[BlackoilPhases::Liquid]) {
121 const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
123 cond[cellIdx].setFreeOil();
127 if (pu.phase_used[BlackoilPhases::Vapour]) {
128 const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
130 cond[cellIdx].setFreeGas();
136 std::vector<std::vector<double>> rho(numPhases);
137 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
138 rho[phaseIdx].resize(numCells);
142 std::vector<double> capPress(numCells*numPhases);
143 std::vector<int> cellIdxArray(numCells);
144 for (
int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
145 cellIdxArray[cellIdx] = cellIdx;
147 props.capPress(numCells, initialState.saturation().data(), cellIdxArray.data(), capPress.data(), NULL);
155 std::vector<std::vector<double> > phasePressure(numPhases);
156 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
157 phasePressure[phaseIdx].resize(numCells);
160 for (
int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
162 assert(pu.phase_used[BlackoilPhases::Liquid]);
164 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
165 phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx];
167 if (pu.phase_used[BlackoilPhases::Aqua]) {
168 const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
169 phasePressure[wpos][cellIdx] =
170 initialState.pressure()[cellIdx]
171 + (capPress[cellIdx*numPhases + opos] - capPress[cellIdx*numPhases + wpos]);
174 if (pu.phase_used[BlackoilPhases::Vapour]) {
175 const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
176 phasePressure[gpos][cellIdx] =
177 initialState.pressure()[cellIdx]
178 + (capPress[cellIdx*numPhases + gpos] - capPress[cellIdx*numPhases + opos]);
183 if (pu.phase_used[BlackoilPhases::Aqua]) {
184 const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
185 const auto& pvtw = props.waterPvt();
186 for (
int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
187 int pvtRegionIdx = pvtRegion[cellIdx];
189 double T = initialState.temperature()[cellIdx];
190 double p = phasePressure[wpos][cellIdx];
191 double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
193 rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
197 if (pu.phase_used[BlackoilPhases::Liquid]) {
198 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
199 const auto& pvto = props.oilPvt();
200 for (
int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
201 int pvtRegionIdx = pvtRegion[cellIdx];
203 double T = initialState.temperature()[cellIdx];
204 double p = phasePressure[opos][cellIdx];
205 double Rs = initialState.gasoilratio()[cellIdx];
206 double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p);
210 b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
213 b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
216 rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
217 if (pu.phase_used[BlackoilPhases::Vapour]) {
218 int gpos = pu.phase_pos[BlackoilPhases::Vapour];
219 rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
224 if (pu.phase_used[BlackoilPhases::Vapour]) {
225 const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
226 const auto& pvtg = props.gasPvt();
227 for (
int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
228 int pvtRegionIdx = pvtRegion[cellIdx];
230 double T = initialState.temperature()[cellIdx];
231 double p = phasePressure[gpos][cellIdx];
232 double Rv = initialState.rv()[cellIdx];
233 double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p);
237 b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
240 b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
242 rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
243 if (pu.phase_used[BlackoilPhases::Liquid]) {
244 int opos = pu.phase_pos[BlackoilPhases::Liquid];
245 rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b;
252 const int num_faces = UgGridHelpers::numFaces(grid);
253 const auto& fc = UgGridHelpers::faceCells(grid);
254 for (
int face = 0; face < num_faces; ++face) {
255 const int c1 = fc(face, 0);
256 const int c2 = fc(face, 1);
257 if (c1 < 0 || c2 < 0) {
261 const int gc1 = (gc == 0) ? c1 : gc[c1];
262 const int gc2 = (gc == 0) ? c2 : gc[c2];
263 const int eq1 = eqlnumData[gc1];
264 const int eq2 = eqlnumData[gc2];
273 const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2));
274 if (maxDp.count(barrierId) == 0) {
275 maxDp[barrierId] = 0.0;
278 for (
int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
279 const double z1 = UgGridHelpers::cellCenterDepth(grid, c1);
280 const double z2 = UgGridHelpers::cellCenterDepth(grid, c2);
282 const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2;
284 const double s1 = initialState.saturation()[numPhases*c1 + phaseIdx];
285 const double s2 = initialState.saturation()[numPhases*c2 + phaseIdx];
287 const double sResid1 = minSat[numPhases*c1 + phaseIdx];
288 const double sResid2 = minSat[numPhases*c2 + phaseIdx];
291 const double p1 = phasePressure[phaseIdx][c1];
292 const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(z1 - z2);
294 if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2))
295 maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
320 template <
class Gr
id>
322 const EclipseState& eclipseState,
324 const std::map<std::pair<int, int>,
double>& maxDp)
326 const SimulationConfig& simulationConfig = eclipseState.getSimulationConfig();
327 std::vector<double> thpres_vals;
328 if (simulationConfig.hasThresholdPressure()) {
329 const ThresholdPressure& thresholdPressure = simulationConfig.getThresholdPressure();
330 const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty(
"EQLNUM");
331 const auto& eqlnumData = eqlnum.getData();
334 const int num_faces = UgGridHelpers::numFaces(grid);
335 const auto& fc = UgGridHelpers::faceCells(grid);
336 const int* gc = UgGridHelpers::globalCell(grid);
337 thpres_vals.resize(num_faces, 0.0);
338 for (
int face = 0; face < num_faces; ++face) {
339 const int c1 = fc(face, 0);
340 const int c2 = fc(face, 1);
341 if (c1 < 0 || c2 < 0) {
345 const int gc1 = (gc == 0) ? c1 : gc[c1];
346 const int gc2 = (gc == 0) ? c2 : gc[c2];
347 const int eq1 = eqlnumData[gc1];
348 const int eq2 = eqlnumData[gc2];
350 if (thresholdPressure.hasRegionBarrier(eq1,eq2)) {
351 if (thresholdPressure.hasThresholdPressure(eq1,eq2)) {
352 thpres_vals[face] = thresholdPressure.getThresholdPressure(eq1,eq2);
358 const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2));
359 if (maxDp.count(barrierId) > 0)
360 thpres_vals[face] = maxDp.at(barrierId);
362 thpres_vals[face] = 0.0;
385 const std::map<std::pair<int, int>,
double>& maxDp)
387 const SimulationConfig& simulationConfig = eclipseState.getSimulationConfig();
388 std::vector<double> thpres_vals;
389 if (simulationConfig.hasThresholdPressure()) {
390 const ThresholdPressure& thresholdPressure = simulationConfig.getThresholdPressure();
391 const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty(
"EQLNUM");
392 const auto& eqlnumData = eqlnum.getData();
396 thpres_vals.resize(nnc.numNNC(), 0.0);
397 for (
size_t i = 0 ; i < nnc.numNNC(); ++i) {
398 const int gc1 = nnc.nncdata()[i].cell1;
399 const int gc2 = nnc.nncdata()[i].cell2;
400 const int eq1 = eqlnumData[gc1];
401 const int eq2 = eqlnumData[gc2];
403 if (thresholdPressure.hasRegionBarrier(eq1,eq2)) {
404 if (thresholdPressure.hasThresholdPressure(eq1,eq2)) {
405 thpres_vals[i] = thresholdPressure.getThresholdPressure(eq1,eq2);
410 const auto barrierId = std::make_pair(eq1, eq2);
411 thpres_vals[i] = maxDp.at(barrierId);
420 #endif // OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
void computeMaxDp(std::map< std::pair< int, int >, double > &maxDp, const Deck &deck, const EclipseState &eclipseState, const Grid &grid, const BlackoilState &initialState, const BlackoilPropertiesFromDeck &props, const double gravity)
Compute the maximum gravity corrected pressure difference of all equilibration regions given a reserv...
Definition: thresholdPressures.hpp:47
std::vector< double > thresholdPressuresNNC(const EclipseState &eclipseState, const NNC &nnc, const std::map< std::pair< int, int >, double > &maxDp)
Get a vector of pressure thresholds from either EclipseState or maxDp (for defaulted values) for all ...
Definition: thresholdPressures.hpp:383
std::vector< double > thresholdPressures(const Deck &, const EclipseState &eclipseState, const Grid &grid, const std::map< std::pair< int, int >, double > &maxDp)
Get a vector of pressure thresholds from EclipseState.
Definition: thresholdPressures.hpp:321