diff --git a/benchmarks/upscale_relperm_benchmark.cpp b/benchmarks/upscale_relperm_benchmark.cpp index e8881a2..54b4441 100644 --- a/benchmarks/upscale_relperm_benchmark.cpp +++ b/benchmarks/upscale_relperm_benchmark.cpp @@ -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 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; diff --git a/examples/upscale_relperm.cpp b/examples/upscale_relperm.cpp index 079df65..859c979 100644 --- a/examples/upscale_relperm.cpp +++ b/examples/upscale_relperm.cpp @@ -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 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; diff --git a/opm/upscaling/RelPermUtils.cpp b/opm/upscaling/RelPermUtils.cpp index ac876da..04c44da 100644 --- a/opm/upscaling/RelPermUtils.cpp +++ b/opm/upscaling/RelPermUtils.cpp @@ -623,4 +623,94 @@ void RelPermUpscaleHelper::calculateMinMaxCapillaryPressure(double dPmin, double } } +void RelPermUpscaleHelper::upscaleCapillaryPressure(std::map& options, + const std::vector& dP) +{ + const double saturationThreshold = atof(options["saturationThreshold"].c_str()); + double largestSaturationInterval = Swor-Swir; + double Ptestvalue = Pcmax; + std::stringstream errstr; + const std::vector& 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 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()); + } +} + } diff --git a/opm/upscaling/RelPermUtils.hpp b/opm/upscaling/RelPermUtils.hpp index 6281504..cb126a8 100644 --- a/opm/upscaling/RelPermUtils.hpp +++ b/opm/upscaling/RelPermUtils.hpp @@ -142,6 +142,13 @@ namespace Opm { //! gravity, linsolver_tolerance void calculateMinMaxCapillaryPressure(double dPmin, double dPmax, std::map& 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& options, + const std::vector& dP); private: //! \brief Perform critical saturation check for a single curve. //! \param[in,out] func Function to check for.