00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
00021 #define OPM_THRESHOLDPRESSURES_HEADER_INCLUDED
00022
00023 #include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
00024 #include <opm/core/props/BlackoilPhases.hpp>
00025
00026 #include <vector>
00027 #include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
00028 #include <opm/parser/eclipse/EclipseState/SimulationConfig/ThresholdPressure.hpp>
00029 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
00030 #include <opm/parser/eclipse/Deck/Deck.hpp>
00031 #include <opm/parser/eclipse/EclipseState/Grid/NNC.hpp>
00032
00033
00034 namespace Opm
00035 {
00046 template <class Grid>
00047 void computeMaxDp(std::map<std::pair<int, int>, double>& maxDp,
00048 const Deck& deck,
00049 const EclipseState& eclipseState,
00050 const Grid& grid,
00051 const BlackoilState& initialState,
00052 const BlackoilPropertiesFromDeck& props,
00053 const double gravity)
00054 {
00055
00056 const PhaseUsage& pu = props.phaseUsage();
00057
00058 const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty("EQLNUM");
00059 const auto& eqlnumData = eqlnum.getData();
00060
00061 const int numPhases = initialState.numPhases();
00062 const int numCells = UgGridHelpers::numCells(grid);
00063 const int numPvtRegions = deck.getKeyword("TABDIMS").getRecord(0).getItem("NTPVT").get< int >(0);
00064
00065
00066 std::vector<double> minSat(numPhases*numCells);
00067 std::vector<double> maxSat(numPhases*numCells);
00068 std::vector<int> allCells(numCells);
00069 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00070 allCells[cellIdx] = cellIdx;
00071 }
00072 props.satRange(numCells, allCells.data(), minSat.data(), maxSat.data());
00073
00074
00075 std::vector<std::vector<double> > surfaceDensity(numPvtRegions);
00076 const auto& densityKw = deck.getKeyword("DENSITY");
00077 for (int regionIdx = 0; regionIdx < numPvtRegions; ++regionIdx) {
00078 surfaceDensity[regionIdx].resize(numPhases);
00079
00080 if (pu.phase_used[BlackoilPhases::Aqua]) {
00081 const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
00082 surfaceDensity[regionIdx][wpos] =
00083 densityKw.getRecord(regionIdx).getItem("WATER").getSIDouble(0);
00084 }
00085
00086 if (pu.phase_used[BlackoilPhases::Liquid]) {
00087 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
00088 surfaceDensity[regionIdx][opos] =
00089 densityKw.getRecord(regionIdx).getItem("OIL").getSIDouble(0);
00090 }
00091
00092 if (pu.phase_used[BlackoilPhases::Vapour]) {
00093 const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
00094 surfaceDensity[regionIdx][gpos] =
00095 densityKw.getRecord(regionIdx).getItem("GAS").getSIDouble(0);
00096 }
00097 }
00098
00099
00100
00101 const int* gc = UgGridHelpers::globalCell(grid);
00102 std::vector<int> pvtRegion(numCells);
00103 const auto& cartPvtRegion = eclipseState.get3DProperties().getIntGridProperty("PVTNUM").getData();
00104 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00105 const int cartCellIdx = gc ? gc[cellIdx] : cellIdx;
00106 pvtRegion[cellIdx] = std::max(0, cartPvtRegion[cartCellIdx] - 1);
00107 }
00108
00109
00110
00111 std::vector<PhasePresence> cond(numCells);
00112 for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) {
00113 if (pu.phase_used[BlackoilPhases::Aqua]) {
00114 const double sw = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Aqua]];
00115 if (sw > 0.0) {
00116 cond[cellIdx].setFreeWater();
00117 }
00118 }
00119
00120 if (pu.phase_used[BlackoilPhases::Liquid]) {
00121 const double so = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Liquid]];
00122 if (so > 0.0) {
00123 cond[cellIdx].setFreeOil();
00124 }
00125 }
00126
00127 if (pu.phase_used[BlackoilPhases::Vapour]) {
00128 const double sg = initialState.saturation()[numPhases*cellIdx + pu.phase_pos[BlackoilPhases::Vapour]];
00129 if (sg > 0.0) {
00130 cond[cellIdx].setFreeGas();
00131 }
00132 }
00133 }
00134
00135
00136 std::vector<std::vector<double>> rho(numPhases);
00137 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00138 rho[phaseIdx].resize(numCells);
00139 }
00140
00141
00142 std::vector<double> capPress(numCells*numPhases);
00143 std::vector<int> cellIdxArray(numCells);
00144 for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
00145 cellIdxArray[cellIdx] = cellIdx;
00146 }
00147 props.capPress(numCells, initialState.saturation().data(), cellIdxArray.data(), capPress.data(), NULL);
00148
00149
00150
00151
00152
00153
00154
00155 std::vector<std::vector<double> > phasePressure(numPhases);
00156 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00157 phasePressure[phaseIdx].resize(numCells);
00158 }
00159
00160 for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
00161
00162 assert(pu.phase_used[BlackoilPhases::Liquid]);
00163
00164 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
00165 phasePressure[opos][cellIdx] = initialState.pressure()[cellIdx];
00166
00167 if (pu.phase_used[BlackoilPhases::Aqua]) {
00168 const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
00169 phasePressure[wpos][cellIdx] =
00170 initialState.pressure()[cellIdx]
00171 + (capPress[cellIdx*numPhases + opos] - capPress[cellIdx*numPhases + wpos]);
00172 }
00173
00174 if (pu.phase_used[BlackoilPhases::Vapour]) {
00175 const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
00176 phasePressure[gpos][cellIdx] =
00177 initialState.pressure()[cellIdx]
00178 + (capPress[cellIdx*numPhases + gpos] - capPress[cellIdx*numPhases + opos]);
00179 }
00180 }
00181
00182
00183 if (pu.phase_used[BlackoilPhases::Aqua]) {
00184 const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
00185 const auto& pvtw = props.waterPvt();
00186 for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
00187 int pvtRegionIdx = pvtRegion[cellIdx];
00188
00189 double T = initialState.temperature()[cellIdx];
00190 double p = phasePressure[wpos][cellIdx];
00191 double b = pvtw.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
00192
00193 rho[wpos][cellIdx] = surfaceDensity[pvtRegionIdx][wpos]*b;
00194 }
00195 }
00196
00197 if (pu.phase_used[BlackoilPhases::Liquid]) {
00198 const int opos = pu.phase_pos[BlackoilPhases::Liquid];
00199 const auto& pvto = props.oilPvt();
00200 for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
00201 int pvtRegionIdx = pvtRegion[cellIdx];
00202
00203 double T = initialState.temperature()[cellIdx];
00204 double p = phasePressure[opos][cellIdx];
00205 double Rs = initialState.gasoilratio()[cellIdx];
00206 double RsSat = pvto.saturatedGasDissolutionFactor(pvtRegionIdx, T, p);
00207
00208 double b;
00209 if (Rs >= RsSat) {
00210 b = pvto.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
00211 }
00212 else {
00213 b = pvto.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rs);
00214 }
00215
00216 rho[opos][cellIdx] = surfaceDensity[pvtRegionIdx][opos]*b;
00217 if (pu.phase_used[BlackoilPhases::Vapour]) {
00218 int gpos = pu.phase_pos[BlackoilPhases::Vapour];
00219 rho[opos][cellIdx] += surfaceDensity[pvtRegionIdx][gpos]*Rs*b;
00220 }
00221 }
00222 }
00223
00224 if (pu.phase_used[BlackoilPhases::Vapour]) {
00225 const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
00226 const auto& pvtg = props.gasPvt();
00227 for (int cellIdx = 0; cellIdx < numCells; ++ cellIdx) {
00228 int pvtRegionIdx = pvtRegion[cellIdx];
00229
00230 double T = initialState.temperature()[cellIdx];
00231 double p = phasePressure[gpos][cellIdx];
00232 double Rv = initialState.rv()[cellIdx];
00233 double RvSat = pvtg.saturatedOilVaporizationFactor(pvtRegionIdx, T, p);
00234
00235 double b;
00236 if (Rv >= RvSat) {
00237 b = pvtg.saturatedInverseFormationVolumeFactor(pvtRegionIdx, T, p);
00238 }
00239 else {
00240 b = pvtg.inverseFormationVolumeFactor(pvtRegionIdx, T, p, Rv);
00241 }
00242 rho[gpos][cellIdx] = surfaceDensity[pvtRegionIdx][gpos]*b;
00243 if (pu.phase_used[BlackoilPhases::Liquid]) {
00244 int opos = pu.phase_pos[BlackoilPhases::Liquid];
00245 rho[gpos][cellIdx] += surfaceDensity[pvtRegionIdx][opos]*Rv*b;
00246 }
00247 }
00248 }
00249
00250
00251
00252 const int num_faces = UgGridHelpers::numFaces(grid);
00253 const auto& fc = UgGridHelpers::faceCells(grid);
00254 for (int face = 0; face < num_faces; ++face) {
00255 const int c1 = fc(face, 0);
00256 const int c2 = fc(face, 1);
00257 if (c1 < 0 || c2 < 0) {
00258
00259 continue;
00260 }
00261 const int gc1 = (gc == 0) ? c1 : gc[c1];
00262 const int gc2 = (gc == 0) ? c2 : gc[c2];
00263 const int eq1 = eqlnumData[gc1];
00264 const int eq2 = eqlnumData[gc2];
00265
00266 if (eq1 == eq2) {
00267
00268 continue;
00269 }
00270
00271
00272
00273 const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2));
00274 if (maxDp.count(barrierId) == 0) {
00275 maxDp[barrierId] = 0.0;
00276 }
00277
00278 for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
00279 const double z1 = UgGridHelpers::cellCenterDepth(grid, c1);
00280 const double z2 = UgGridHelpers::cellCenterDepth(grid, c2);
00281
00282 const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2;
00283
00284 const double s1 = initialState.saturation()[numPhases*c1 + phaseIdx];
00285 const double s2 = initialState.saturation()[numPhases*c2 + phaseIdx];
00286
00287 const double sResid1 = minSat[numPhases*c1 + phaseIdx];
00288 const double sResid2 = minSat[numPhases*c2 + phaseIdx];
00289
00290
00291 const double p1 = phasePressure[phaseIdx][c1];
00292 const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(z1 - z2);
00293
00294 if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2))
00295 maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));
00296 }
00297 }
00298 }
00299
00317
00318
00319
00320 template <class Grid>
00321 std::vector<double> thresholdPressures(const Deck& ,
00322 const EclipseState& eclipseState,
00323 const Grid& grid,
00324 const std::map<std::pair<int, int>, double>& maxDp)
00325 {
00326 const SimulationConfig& simulationConfig = eclipseState.getSimulationConfig();
00327 std::vector<double> thpres_vals;
00328 if (simulationConfig.hasThresholdPressure()) {
00329 const ThresholdPressure& thresholdPressure = simulationConfig.getThresholdPressure();
00330 const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty("EQLNUM");
00331 const auto& eqlnumData = eqlnum.getData();
00332
00333
00334 const int num_faces = UgGridHelpers::numFaces(grid);
00335 const auto& fc = UgGridHelpers::faceCells(grid);
00336 const int* gc = UgGridHelpers::globalCell(grid);
00337 thpres_vals.resize(num_faces, 0.0);
00338 for (int face = 0; face < num_faces; ++face) {
00339 const int c1 = fc(face, 0);
00340 const int c2 = fc(face, 1);
00341 if (c1 < 0 || c2 < 0) {
00342
00343 continue;
00344 }
00345 const int gc1 = (gc == 0) ? c1 : gc[c1];
00346 const int gc2 = (gc == 0) ? c2 : gc[c2];
00347 const int eq1 = eqlnumData[gc1];
00348 const int eq2 = eqlnumData[gc2];
00349
00350 if (thresholdPressure.hasRegionBarrier(eq1,eq2)) {
00351 if (thresholdPressure.hasThresholdPressure(eq1,eq2)) {
00352 thpres_vals[face] = thresholdPressure.getThresholdPressure(eq1,eq2);
00353 }
00354 else {
00355
00356
00357
00358 const auto barrierId = std::make_pair(std::min(eq1, eq2), std::max(eq1, eq2));
00359 if (maxDp.count(barrierId) > 0)
00360 thpres_vals[face] = maxDp.at(barrierId);
00361 else
00362 thpres_vals[face] = 0.0;
00363 }
00364 }
00365
00366 }
00367 }
00368 return thpres_vals;
00369 }
00370
00383 std::vector<double> thresholdPressuresNNC(const EclipseState& eclipseState,
00384 const NNC& nnc,
00385 const std::map<std::pair<int, int>, double>& maxDp)
00386 {
00387 const SimulationConfig& simulationConfig = eclipseState.getSimulationConfig();
00388 std::vector<double> thpres_vals;
00389 if (simulationConfig.hasThresholdPressure()) {
00390 const ThresholdPressure& thresholdPressure = simulationConfig.getThresholdPressure();
00391 const auto& eqlnum = eclipseState.get3DProperties().getIntGridProperty("EQLNUM");
00392 const auto& eqlnumData = eqlnum.getData();
00393
00394
00395
00396 thpres_vals.resize(nnc.numNNC(), 0.0);
00397 for (size_t i = 0 ; i < nnc.numNNC(); ++i) {
00398 const int gc1 = nnc.nncdata()[i].cell1;
00399 const int gc2 = nnc.nncdata()[i].cell2;
00400 const int eq1 = eqlnumData[gc1];
00401 const int eq2 = eqlnumData[gc2];
00402
00403 if (thresholdPressure.hasRegionBarrier(eq1,eq2)) {
00404 if (thresholdPressure.hasThresholdPressure(eq1,eq2)) {
00405 thpres_vals[i] = thresholdPressure.getThresholdPressure(eq1,eq2);
00406 } else {
00407
00408
00409
00410 const auto barrierId = std::make_pair(eq1, eq2);
00411 thpres_vals[i] = maxDp.at(barrierId);
00412 }
00413 }
00414 }
00415 }
00416 return thpres_vals;
00417 }
00418 }
00419
00420 #endif // OPM_THRESHOLDPRESSURES_HEADER_INCLUDED