changed: put upscaling of capillary pressure in a function

- 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 ddf8819f94
commit 0f87352e92
4 changed files with 101 additions and 181 deletions
+2 -91
View File
@@ -452,100 +452,11 @@ try
* will be used afterwards for accessing the tabulated values.
*/
helper.upscaleCapillaryPressure(options, dP);
double largestSaturationInterval = helper.Swor-helper.Swir;
double Ptestvalue = helper.Pcmax;
while (largestSaturationInterval > (helper.Swor-helper.Swir)/500.0) {
// cout << Ptestvalue << endl;
if (helper.Pcmax == helper.Pcmin) {
// This is a dummy situation, we go through once and then
// we are finished (this will be triggered by zero permeability)
Ptestvalue = helper.Pcmin;
largestSaturationInterval = 0;
}
else if (helper.WaterSaturationVsCapPressure.getSize() == 0) {
/* No data values previously computed */
Ptestvalue = helper.Pcmax;
}
else if (helper.WaterSaturationVsCapPressure.getSize() == 1) {
/* If only one point has been computed, it was for Pcmax. So now
do Pcmin */
Ptestvalue = helper.Pcmin;
}
else {
/* Search for largest saturation interval in which there are no
computed saturation points (and estimate the capillary pressure
that will fall in the center of this saturation interval)
*/
pair<double,double> SatDiff = helper.WaterSaturationVsCapPressure.getMissingX();
Ptestvalue = SatDiff.first;
largestSaturationInterval = SatDiff.second;
}
// Check for saneness of Ptestvalue:
if (std::isnan(Ptestvalue) || std::isinf(Ptestvalue)) {
if (helper.isMaster) cerr << "ERROR: Ptestvalue was inf or nan" << endl;
break; // Jump out of while-loop, just print out the results
// up to now and exit the program
}
double waterVolume = 0.0;
for (unsigned int i = 0; i < ecl_idx.size(); ++i) {
unsigned int cell_idx = ecl_idx[i];
double WaterSaturationCell = 0.0;
if (helper.satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero
double PtestvalueCell;
if (includeGravity) {
PtestvalueCell = Ptestvalue - dP[cell_idx];
}
else {
PtestvalueCell = Ptestvalue;
}
if (! helper.anisotropic_input ) {
double Jvalue = sqrt(helper.perms[0][cell_idx] * milliDarcyToSqMetre /helper.poros[cell_idx]) * PtestvalueCell;
//cout << "JvalueCell: " << Jvalue << endl;
WaterSaturationCell
= helper.InvJfunctions[int(helper.satnums[cell_idx])-1].evaluate(Jvalue);
}
else { // anisotropic_input, then we do not do J-function-scaling
WaterSaturationCell = helper.SwPcfunctions[int(helper.satnums[cell_idx])-1].evaluate(PtestvalueCell);
//cout << Ptestvalue << "\t" << helper.WaterSaturationCell << endl;
}
}
waterVolume += WaterSaturationCell * helper.cellPoreVolumes[cell_idx];
}
helper.WaterSaturationVsCapPressure.addPair(Ptestvalue, waterVolume/helper.poreVolume);
}
// cout << WaterSaturationVsCapPressure.toString();
// Now, it may happen that we have a large number of cells, and
// some cells with near zero poro and perm. This may cause that
// Pcmax has been estimated so high that it does not affect Sw
// within machine precision, and then we need to truncate the
// largest Pc values:
helper.WaterSaturationVsCapPressure.chopFlatEndpoints(saturationThreshold);
// Now we can also invert the upscaled water saturation
// (it should be monotonic)
if (!helper.WaterSaturationVsCapPressure.isStrictlyMonotone()) {
if (helper.isMaster) {
cerr << "Error: Upscaled water saturation not strictly monotone in capillary pressure." << endl;
cerr << " Unphysical input data, exiting." << endl;
cerr << " Trying to dump " << helper.saturationstring << " vs Pc to file swvspc_debug.txt for inspection" << endl;
ofstream outfile;
outfile.open("swvspc_debug.txt", ios::out | ios::trunc);
outfile << "# Pc " << helper.saturationstring << endl;
outfile << helper.WaterSaturationVsCapPressure.toString();
outfile.close();
}
usageandexit();
}
MonotCubicInterpolator CapPressureVsWaterSaturation(helper.WaterSaturationVsCapPressure.get_fVector(),
helper.WaterSaturationVsCapPressure.get_xVector());
clock_t start_upscaling = clock();
/*****************************************************************************
@@ -632,7 +543,7 @@ try
// Should "I" (mpi-wise) compute this pressure point?
if (helper.node_vs_pressurepoint[pointidx] == mpi_rank) {
Ptestvalue = helper.pressurePoints[pointidx];
double Ptestvalue = helper.pressurePoints[pointidx];
double accPhasePerm = 0.0;
double accPhase2Perm = 0.0;
+2 -90
View File
@@ -672,95 +672,7 @@ try
* will be used afterwards for accessing the tabulated values.
*/
double largestSaturationInterval = helper.Swor-helper.Swir;
double Ptestvalue = helper.Pcmax;
while (largestSaturationInterval > (helper.Swor-helper.Swir)/500.0) {
// cout << Ptestvalue << endl;
if (helper.Pcmax == helper.Pcmin) {
// This is a dummy situation, we go through once and then
// we are finished (this will be triggered by zero permeability)
Ptestvalue = helper.Pcmin;
largestSaturationInterval = 0;
}
else if (helper.WaterSaturationVsCapPressure.getSize() == 0) {
/* No data values previously computed */
Ptestvalue = helper.Pcmax;
}
else if (helper.WaterSaturationVsCapPressure.getSize() == 1) {
/* If only one point has been computed, it was for Pcmax. So now
do Pcmin */
Ptestvalue = helper.Pcmin;
}
else {
/* Search for largest saturation interval in which there are no
computed saturation points (and estimate the capillary pressure
that will fall in the center of this saturation interval)
*/
pair<double,double> SatDiff = helper.WaterSaturationVsCapPressure.getMissingX();
Ptestvalue = SatDiff.first;
largestSaturationInterval = SatDiff.second;
}
// Check for saneness of Ptestvalue:
if (std::isnan(Ptestvalue) || std::isinf(Ptestvalue)) {
if (helper.isMaster) cerr << "ERROR: Ptestvalue was inf or nan" << endl;
break; // Jump out of while-loop, just print out the results
// up to now and exit the program
}
double waterVolume = 0.0;
for (unsigned int i = 0; i < ecl_idx.size(); ++i) {
unsigned int cell_idx = ecl_idx[i];
double waterSaturationCell = 0.0;
if (helper.satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero
double PtestvalueCell;
if (includeGravity) {
PtestvalueCell = Ptestvalue - dP[cell_idx];
}
else {
PtestvalueCell = Ptestvalue;
}
if (! helper.anisotropic_input ) {
double Jvalue = sqrt(helper.perms[0][cell_idx] * milliDarcyToSqMetre /helper.poros[cell_idx]) * PtestvalueCell;
//cout << "JvalueCell: " << Jvalue << endl;
waterSaturationCell
= helper.InvJfunctions[int(helper.satnums[cell_idx])-1].evaluate(Jvalue);
}
else { // anisotropic_input, then we do not do J-function-scaling
waterSaturationCell = helper.SwPcfunctions[int(helper.satnums[cell_idx])-1].evaluate(PtestvalueCell);
//cout << Ptestvalue << "\t" << waterSaturationCell << endl;
}
}
waterVolume += waterSaturationCell * helper.cellPoreVolumes[cell_idx];
}
helper.WaterSaturationVsCapPressure.addPair(Ptestvalue, waterVolume/helper.poreVolume);
}
// cout << WaterSaturationVsCapPressure.toString();
// Now, it may happen that we have a large number of cells, and
// some cells with near zero poro and perm. This may cause that
// Pcmax has been estimated so high that it does not affect Sw
// within machine precision, and then we need to truncate the
// largest Pc values:
helper.WaterSaturationVsCapPressure.chopFlatEndpoints(saturationThreshold);
// Now we can also invert the upscaled water saturation
// (it should be monotonic)
if (!helper.WaterSaturationVsCapPressure.isStrictlyMonotone()) {
if (helper.isMaster) {
cerr << "Error: Upscaled water saturation not strictly monotone in capillary pressure." << endl;
cerr << " Unphysical input data, exiting." << endl;
cerr << " Trying to dump " << helper.saturationstring << " vs Pc to file swvspc_debug.txt for inspection" << endl;
ofstream outfile;
outfile.open("swvspc_debug.txt", ios::out | ios::trunc);
outfile << "# Pc " << helper.saturationstring << endl;
outfile << helper.WaterSaturationVsCapPressure.toString();
outfile.close();
}
usageandexit();
}
helper.upscaleCapillaryPressure(options, dP);
MonotCubicInterpolator CapPressureVsWaterSaturation(helper.WaterSaturationVsCapPressure.get_fVector(),
helper.WaterSaturationVsCapPressure.get_xVector());
@@ -848,7 +760,7 @@ try
// Should "I" (mpi-wise) compute this pressure point?
if (helper.node_vs_pressurepoint[pointidx] == mpi_rank) {
Ptestvalue = helper.pressurePoints[pointidx];
double Ptestvalue = helper.pressurePoints[pointidx];
double accPhasePerm = 0.0;
double accPhase2Perm = 0.0;
+90
View File
@@ -623,4 +623,94 @@ void RelPermUpscaleHelper::calculateMinMaxCapillaryPressure(double dPmin, double
}
}
void RelPermUpscaleHelper::upscaleCapillaryPressure(std::map<std::string,std::string>& options,
const std::vector<double>& dP)
{
const double saturationThreshold = atof(options["saturationThreshold"].c_str());
double largestSaturationInterval = Swor-Swir;
double Ptestvalue = Pcmax;
std::stringstream errstr;
const std::vector<int>& ecl_idx = upscaler.grid().globalCell();
const double milliDarcyToSqMetre =
Opm::unit::convert::to(1.0*Opm::prefix::milli*Opm::unit::darcy,
Opm::unit::square(Opm::unit::meter));
while (largestSaturationInterval > (Swor-Swir)/500.0) {
if (Pcmax == Pcmin) {
// This is a dummy situation, we go through once and then
// we are finished (this will be triggered by zero permeability)
Ptestvalue = Pcmin;
largestSaturationInterval = 0;
}
else if (WaterSaturationVsCapPressure.getSize() == 0) {
/* No data values previously computed */
Ptestvalue = Pcmax;
}
else if (WaterSaturationVsCapPressure.getSize() == 1) {
/* If only one point has been computed, it was for Pcmax. So now
do Pcmin */
Ptestvalue = Pcmin;
}
else {
/* Search for largest saturation interval in which there are no
computed saturation points (and estimate the capillary pressure
that will fall in the center of this saturation interval)
*/
std::pair<double,double> SatDiff = WaterSaturationVsCapPressure.getMissingX();
Ptestvalue = SatDiff.first;
largestSaturationInterval = SatDiff.second;
}
// Check for saneness of Ptestvalue:
if (std::isnan(Ptestvalue) || std::isinf(Ptestvalue)) {
errstr << "ERROR: Ptestvalue was inf or nan" << std::endl;
break; // Jump out of while-loop, just print out the results
// up to now and exit the program
}
double waterVolume = 0.0;
for (size_t i = 0; i < ecl_idx.size(); ++i) {
unsigned int cell_idx = ecl_idx[i];
double waterSaturationCell = 0.0;
if (satnums[cell_idx] > 0) { // handle "no rock" cells with satnum zero
double PtestvalueCell = Ptestvalue;
if (!dP.empty())
PtestvalueCell -= dP[cell_idx];
if (!anisotropic_input) {
double Jvalue = sqrt(perms[0][cell_idx] * milliDarcyToSqMetre / poros[cell_idx]) * PtestvalueCell;
waterSaturationCell = InvJfunctions[int(satnums[cell_idx])-1].evaluate(Jvalue);
}
else // anisotropic_input, then we do not do J-function-scaling
waterSaturationCell = SwPcfunctions[int(satnums[cell_idx])-1].evaluate(PtestvalueCell);
}
waterVolume += waterSaturationCell * cellPoreVolumes[cell_idx];
}
WaterSaturationVsCapPressure.addPair(Ptestvalue, waterVolume/poreVolume);
}
// Now, it may happen that we have a large number of cells, and
// some cells with near zero poro and perm. This may cause that
// Pcmax has been estimated so high that it does not affect Sw
// within machine precision, and then we need to truncate the
// largest Pc values:
WaterSaturationVsCapPressure.chopFlatEndpoints(saturationThreshold);
// Now we can also invert the upscaled water saturation
// (it should be monotonic)
if (!WaterSaturationVsCapPressure.isStrictlyMonotone()) {
errstr << "Error: Upscaled water saturation not strictly monotone in capillary pressure." << std::endl
<< " Unphysical input data, exiting." << std::endl
<< " Trying to dump " << saturationstring << " vs Pc to file swvspc_debug.txt for inspection";
if (isMaster) {
std::ofstream outfile;
outfile.open("swvspc_debug.txt", std::ios::out | std::ios::trunc);
outfile << "# Pc " << saturationstring << std::endl;
outfile << WaterSaturationVsCapPressure.toString();
outfile.close();
}
throw std::runtime_error(errstr.str());
}
}
}
+7
View File
@@ -142,6 +142,13 @@ namespace Opm {
//! gravity, linsolver_tolerance
void calculateMinMaxCapillaryPressure(double dPmin, double dPmax,
std::map<std::string,std::string>& options);
//! \brief Upscale capillary pressure.
//! \param[in] Options structure.
//! \param[in] dP Vector with cell pressure gradients. Pass empty vector for no gravity effects.
//! \details Uses the following options: saturationThreshold
void upscaleCapillaryPressure(std::map<std::string,std::string>& options,
const std::vector<double>& dP);
private:
//! \brief Perform critical saturation check for a single curve.
//! \param[in,out] func Function to check for.