OilfieldToolsNet/opm-upscaling
4
0
Fork 0
mirror of https://github.com/OPM/opm-upscaling.git synced 2026-08-26 21:17:11 -05:00
Code Issues Packages Projects Releases Wiki Activity

changed: move calculation of cell center pressure gradients to a function

Browse Source 
- was duplicated across upscale_relperm / upscale_relperm_benchmark
This commit is contained in:
Arne Morten Kvarving
2015-09-14 11:34:31 +02:00
parent 15d47ef70e
commit 95a6461107
4 changed files with 75 additions and 114 deletions
Download Patch File Download Diff File Expand all files Collapse all files
benchmarks/upscale_relperm_benchmark.cpp
+3 -57
View File
@@ -372,8 +372,6 @@ try
// Input for surfaceTension is dynes/cm
// SI units are Joules/square metre
const double surfaceTension = atof(options["surfaceTension"].c_str()) * 1e-3; // multiply with 10^-3 to obtain SI units
const double waterDensity = atof(options["waterDensity"].c_str());
const double oilDensity = atof(options["oilDensity"].c_str());
const bool includeGravity = (fabs(gravity) > DBL_MIN); // true for non-zero gravity
//const int outputprecision = atoi(options["outputprecision"].c_str());
@@ -446,63 +444,11 @@ try
/* If gravity is to be included, calculate z-values of every cell: */
if (includeGravity) {
// height of model is calculated as the average of the z-values at the top layer
// This calculation makes assumption on the indexing of cells in the grid, going from bottom to top.
double modelHeight = 0;
for (unsigned int zIdx = (4 * res[0] * res[1] * (2*res[2]-1)); zIdx < helper.zcorns.size(); ++zIdx) {
modelHeight += helper.zcorns[zIdx] / (4*res[0]*res[1]);
}
// We assume that the spatial units in the grid file is in centimetres,
// so we divide by 100 to get to metres.
modelHeight = modelHeight/100.0;
// Input water and oil density is given in g/cm3, we convert it to kg/m3 (SI)
// by multiplying with 1000.
double dRho = (waterDensity-oilDensity) * 1000; // SI unit (kg/m3)
// Calculating difference in capillary pressure for all cells
dP = vector<double>(helper.satnums.size(), 0);
for (unsigned int cellIdx = 0; cellIdx < helper.satnums.size(); ++cellIdx) {
int i,j,k; // Position of cell in cell hierarchy
vector<int> zIndices(8,0); // 8 corners with 8 heights
int horIdx = (cellIdx+1) - int(std::floor(((double)(cellIdx+1))/((double)(res[0]*res[1]))))*res[0]*res[1]; // index in the corresponding horizon
if (horIdx == 0) {
horIdx = res[0]*res[1];
}
i = horIdx - int(std::floor(((double)horIdx)/((double)res[0])))*res[0];
if (i == 0) {
i = res[0];
}
j = (horIdx-i)/res[0]+1;
k = ((cellIdx+1)-res[0]*(j-1)-1)/(res[0]*res[1])+1;
int zBegin = 8*res[0]*res[1]*(k-1); // indices of Z-values of bottom
int level2 = 4*res[0]*res[1]; // number of z-values in one horizon
zIndices[0] = zBegin + 4*res[0]*(j-1)+2*i-1;
zIndices[1] = zBegin + 4*res[0]*(j-1)+2*i;
zIndices[2] = zBegin + 2*res[0]*(2*j-1)+2*i;
zIndices[3] = zBegin + 2*res[0]*(2*j-1)+2*i-1;
zIndices[4] = zBegin + level2 + 4*res[0]*(j-1)+2*i-1;
zIndices[5] = zBegin + level2 + 4*res[0]*(j-1)+2*i;
zIndices[6] = zBegin + level2 + 2*res[0]*(2*j-1)+2*i;
zIndices[7] = zBegin + level2 + 2*res[0]*(2*j-1)+2*i-1;
double cellDepth = 0;
for (unsigned int corner = 0; corner < 8; ++corner) {
cellDepth += helper.zcorns[zIndices[corner]-1] / 8.0;
}
// cellDepth is in cm, convert to m by dividing by 100
cellDepth = cellDepth / 100.0;
dP[cellIdx] = dRho * gravity * (cellDepth-modelHeight/2.0);
// assume distances in grid are given in cm.
dPmin = min(dPmin,dP[cellIdx]);
dPmax = max(dPmax,dP[cellIdx]);
}
dP = helper.calculateCellPressureGradients(res, options);
dPmin = *std::min_element(dP.begin(), dP.end());
dPmax = *std::max_element(dP.begin(), dP.end());
}
/******************************************************************************
* Step 5:
* Go through each cell and calculate the minimum and