All Classes Namespaces Files Functions Variables Typedefs Enumerator Pages
thresholdPressures.hpp
1 /*
2  Copyright 2014 SINTEF ICT, Applied Mathematics.
3 
4  This file is part of the Open Porous Media project (OPM).
5 
6  OPM is free software: you can redistribute it and/or modify
7  it under the terms of the GNU General Public License as published by
8  the Free Software Foundation, either version 3 of the License, or
9  (at your option) any later version.
10 
11  OPM is distributed in the hope that it will be useful,
12  but WITHOUT ANY WARRANTY; without even the implied warranty of
13  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  GNU General Public License for more details.
15 
16  You should have received a copy of the GNU General Public License
17  along with OPM. If not, see <http://www.gnu.org/licenses/>.
18 */
19 
20 #ifndef OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
21 #define OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
22 
23 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
24 #include <opm/core/props/BlackoilPhases.hpp>
25 
26 #include <vector>
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>
32 
33 
34 namespace Opm
35 {
46 template <class Grid>
47 void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
48  const Deck& deck,
49  const EclipseState& eclipseState,
50  const Grid& grid,
51  const BlackoilState& initialState,
52  const BlackoilPropertiesFromDeck& props,
53  const double gravity)
54 {
55 
56  const PhaseUsage& pu = props.phaseUsage();
57 
58  const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty("EQLNUM");
59  const auto& eqlnumData = eqlnum.getData();
60 
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);
64 
65  // retrieve the minimum (residual!?) and the maximum saturations for all cells
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;
71  }
72  props.satRange(numCells, allCells.data(), minSat.data(), maxSat.data());
73 
74  // retrieve the surface densities
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);
79 
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);
84  }
85 
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);
90  }
91 
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);
96  }
97  }
98 
99  // retrieve the PVT region of each cell. note that we need c++ instead of
100  // Fortran indices.
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);
107  }
108 
109  // compute the initial "phase presence" of each cell (required to calculate
110  // the inverse formation volume factors
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]];
115  if (sw > 0.0) {
116  cond[cellIdx].setFreeWater();
117  }
118  }
119 
120  if (pu.phase_used[BlackoilPhases::Liquid]) {
121  const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
122  if (so > 0.0) {
123  cond[cellIdx].setFreeOil();
124  }
125  }
126 
127  if (pu.phase_used[BlackoilPhases::Vapour]) {
128  const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
129  if (sg > 0.0) {
130  cond[cellIdx].setFreeGas();
131  }
132  }
133  }
134 
135  // calculate the initial fluid densities for the gravity correction.
136  std::vector<std::vector<double>> rho(numPhases);
137  for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
138  rho[phaseIdx].resize(numCells);
139  }
140 
141  // compute the capillary pressures of the active phases
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;
146  }
147  props.capPress(numCells, initialState.saturation().data(), cellIdxArray.data(), capPress.data(), NULL);
148 
149  // compute the absolute pressure of each active phase: for some reason, E100
150  // defines the capillary pressure for the water phase as p_o - p_w while it
151  // uses p_g - p_o for the gas phase. (it would be more consistent to use the
152  // oil pressure as reference for both the other phases.) probably this is
153  // done to always have a positive number for the capillary pressure (as long
154  // as the medium is hydrophilic)
155  std::vector<std::vector<double> > phasePressure(numPhases);
156  for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
157  phasePressure[phaseIdx].resize(numCells);
158  }
159 
160  for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
161  // we currently hard-code the oil phase as the reference phase!
162  assert(pu.phase_used[BlackoilPhases::Liquid]);
163 
164  const int opos = pu.phase_pos[BlackoilPhases::Liquid];
165  phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx];
166 
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]);
172  }
173 
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]);
179  }
180  }
181 
182  // calculate the densities of the active phases for each cell
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];
188 
189  double T = initialState.temperature()[cellIdx];
190  double p = phasePressure[wpos][cellIdx];
191  double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
192 
193  rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
194  }
195  }
196 
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];
202 
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);
207 
208  double b;
209  if (Rs >= RsSat) {
210  b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
211  }
212  else {
213  b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
214  }
215 
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;
220  }
221  }
222  }
223 
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];
229 
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);
234 
235  double b;
236  if (Rv >= RvSat) {
237  b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
238  }
239  else {
240  b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
241  }
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;
246  }
247  }
248  }
249 
250  // Calculate the maximum pressure potential difference between all PVT region
251  // transitions of the initial solution.
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) {
258  // Boundary face, skip this.
259  continue;
260  }
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];
265 
266  if (eq1 == eq2) {
267  // not an equilibration region boundary. skip this.
268  continue;
269  }
270 
271  // update the maximum pressure potential difference between the two
272  // regions
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;
276  }
277 
278  for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
279  const double z1 = UgGridHelpers::cellCenterDepth(grid, c1);
280  const double z2 = UgGridHelpers::cellCenterDepth(grid, c2);
281 
282  const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2;
283 
284  const double s1 = initialState.saturation()[numPhases*c1 + phaseIdx];
285  const double s2 = initialState.saturation()[numPhases*c2 + phaseIdx];
286 
287  const double sResid1 = minSat[numPhases*c1 + phaseIdx];
288  const double sResid2 = minSat[numPhases*c2 + phaseIdx];
289 
290  // compute gravity corrected pressure potentials at the average depth
291  const double p1 = phasePressure[phaseIdx][c1];
292  const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(z1 - z2);
293 
294  if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2))
295  maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
296  }
297  }
298 }
299 
317 
318 
319 
320  template <class Grid>
321  std::vector<double> thresholdPressures(const Deck& /* deck */,
322  const EclipseState& eclipseState,
323  const Grid& grid,
324  const std::map<std::pair<int, int>, double>& maxDp)
325  {
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();
332 
333  // Set threshold pressure values for each cell face.
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) {
342  // Boundary face, skip it.
343  continue;
344  }
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];
349 
350  if (thresholdPressure.hasRegionBarrier(eq1,eq2)) {
351  if (thresholdPressure.hasThresholdPressure(eq1,eq2)) {
352  thpres_vals[face] = thresholdPressure.getThresholdPressure(eq1,eq2);
353  }
354  else {
355  // set the threshold pressure for faces of PVT regions where the third item
356  // has been defaulted to the maximum pressure potential difference between
357  // these regions
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);
361  else
362  thpres_vals[face] = 0.0;
363  }
364  }
365 
366  }
367  }
368  return thpres_vals;
369  }
370 
383  std::vector<double> thresholdPressuresNNC(const EclipseState& eclipseState,
384  const NNC& nnc,
385  const std::map<std::pair<int, int>, double>& maxDp)
386  {
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();
393 
394  // Set values for each NNC
395 
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];
402 
403  if (thresholdPressure.hasRegionBarrier(eq1,eq2)) {
404  if (thresholdPressure.hasThresholdPressure(eq1,eq2)) {
405  thpres_vals[i] = thresholdPressure.getThresholdPressure(eq1,eq2);
406  } else {
407  // set the threshold pressure for NNC of PVT regions where the third item
408  // has been defaulted to the maximum pressure potential difference between
409  // these regions
410  const auto barrierId = std::make_pair(eq1, eq2);
411  thpres_vals[i] = maxDp.at(barrierId);
412  }
413  }
414  }
415  }
416  return thpres_vals;
417  }
418 }
419 
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