From f1f1df5c4f1ebc36a31185207f7ecef00dced28a Mon Sep 17 00:00:00 2001 From: Arne Morten Kvarving Date: Wed, 1 Jul 2015 16:10:47 +0200 Subject: [PATCH] refactor: move phase perm to relperm calculations to helper class - avoids duplication per phase - avoids duplication in upscale_relperm / upscale_relperm_benchmark --- benchmarks/upscale_relperm_benchmark.cpp | 88 +---------------- examples/upscale_relperm.cpp | 115 ++++------------------- opm/upscaling/RelPermUtils.cpp | 44 +++++++++ opm/upscaling/RelPermUtils.hpp | 6 ++ 4 files changed, 72 insertions(+), 181 deletions(-) diff --git a/benchmarks/upscale_relperm_benchmark.cpp b/benchmarks/upscale_relperm_benchmark.cpp index 3dc2ee0..aef96fa 100644 --- a/benchmarks/upscale_relperm_benchmark.cpp +++ b/benchmarks/upscale_relperm_benchmark.cpp @@ -1394,90 +1394,10 @@ try * Step 8c: Make relperm values from phaseperms * (only master node can do this) */ - - vector > RelPermValues; // voigtIdx is first index. - for (int voigtIdx=0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - vector tmp; - RelPermValues.push_back(tmp); - } - if (helper.isMaster) { - // Loop over all pressure points - for (int idx=0; idx < helper.points; ++idx) { - Matrix phasePermTensor = zeroMatrix; - zero(phasePermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - setVoigtValue(phasePermTensor, voigtIdx, helper.PhasePerm[0][idx][voigtIdx]); - } - //cout << phasePermTensor << endl; - Matrix relPermTensor = zeroMatrix; - // relPermTensor = phasePermTensor; - // relPermTensor *= permTensorInv; - prod(phasePermTensor, helper.permTensorInv, relPermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - RelPermValues[voigtIdx].push_back(getVoigtValue(relPermTensor, voigtIdx)); - } - //cout << relPermTensor << endl; - } - } - - vector > RelPermValues2; // voigtIdx is first index. - if (helper.upscaleBothPhases) { - for (int voigtIdx=0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - vector tmp; - RelPermValues2.push_back(tmp); - } - if (helper.isMaster) { - // Loop over all pressure points - for (int idx=0; idx < helper.points; ++idx) { - Matrix phasePermTensor = zeroMatrix; - zero(phasePermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - setVoigtValue(phasePermTensor, voigtIdx, helper.PhasePerm[1][idx][voigtIdx]); - } - //cout << phasePermTensor << endl; - Matrix relPermTensor = zeroMatrix; - // relPermTensor = phasePermTensor; - // relPermTensor *= permTensorInv; - prod(phasePermTensor, helper.permTensorInv, relPermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - RelPermValues2[voigtIdx].push_back(getVoigtValue(relPermTensor, voigtIdx)); - } - //cout << relPermTensor << endl; - } - } - } - - // If doEclipseCheck, critical saturation points should be specified by 0 relperm - // Numerical errors and maxpermcontrast violate this even if the input has specified - // these points - if (helper.isMaster) { - if (helper.doEclipseCheck) { - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - int minidx; - if (RelPermValues[voigtIdx][0] < RelPermValues[voigtIdx][helper.points-1]) minidx = 0; else minidx = helper.points-1; - if (RelPermValues[voigtIdx][minidx] < helper.critRelpThresh) { - RelPermValues[voigtIdx][minidx] = 0.0; - } - else { - cerr << "Minimum upscaled relperm value is " << RelPermValues[voigtIdx][minidx] << ", larger than critRelpermThresh." << endl - << "(voigtidx = " << voigtIdx << ")" << endl; - usageandexit(); - } - if (helper.upscaleBothPhases) { - if (RelPermValues2[voigtIdx][0] < RelPermValues2[voigtIdx][helper.points-1]) minidx = 0; else minidx = helper.points-1; - if (RelPermValues2[voigtIdx][minidx] < helper.critRelpThresh) { - RelPermValues2[voigtIdx][minidx] = 0.0; - } - else { - cerr << "Minimum upscaled relperm value for phase 2 is " << RelPermValues2[voigtIdx][minidx] << endl - << ", larger than critRelpermThresh.(voigtidx = " << voigtIdx << ")" << endl; - usageandexit(); - } - } - } - } - } - + std::array>,2> RelPermValues; + RelPermValues[0] = helper.getRelPerm(0); + if (helper.upscaleBothPhases) + RelPermValues[1] = helper.getRelPerm(1); /********************************************************************************* * Step 9 - Benchmark version diff --git a/examples/upscale_relperm.cpp b/examples/upscale_relperm.cpp index 2f6a502..53e04ba 100644 --- a/examples/upscale_relperm.cpp +++ b/examples/upscale_relperm.cpp @@ -1532,89 +1532,10 @@ try * Step 8c: Make relperm values from phaseperms * (only master node can do this) */ - - vector > RelPermValues; // voigtIdx is first index. - for (int voigtIdx=0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - vector tmp; - RelPermValues.push_back(tmp); - } - if (helper.isMaster) { - // Loop over all pressure points - for (int idx=0; idx < helper.points; ++idx) { - Matrix phasePermTensor = zeroMatrix; - zero(phasePermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - setVoigtValue(phasePermTensor, voigtIdx, helper.PhasePerm[0][idx][voigtIdx]); - } - //cout << phasePermTensor << endl; - Matrix relPermTensor = zeroMatrix; - // relPermTensor = phasePermTensor; - // relPermTensor *= permTensorInv; - prod(phasePermTensor, helper.permTensorInv, relPermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - RelPermValues[voigtIdx].push_back(getVoigtValue(relPermTensor, voigtIdx)); - } - //cout << relPermTensor << endl; - } - } - - vector > RelPermValues2; // voigtIdx is first index. - if (helper.upscaleBothPhases) { - for (int voigtIdx=0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - vector tmp; - RelPermValues2.push_back(tmp); - } - if (helper.isMaster) { - // Loop over all pressure points - for (int idx=0; idx < helper.points; ++idx) { - Matrix phasePermTensor = zeroMatrix; - zero(phasePermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - setVoigtValue(phasePermTensor, voigtIdx, helper.PhasePerm[1][idx][voigtIdx]); - } - //cout << phasePermTensor << endl; - Matrix relPermTensor = zeroMatrix; - // relPermTensor = phasePermTensor; - // relPermTensor *= permTensorInv; - prod(phasePermTensor, helper.permTensorInv, relPermTensor); - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - RelPermValues2[voigtIdx].push_back(getVoigtValue(relPermTensor, voigtIdx)); - } - //cout << relPermTensor << endl; - } - } - } - - // If doEclipseCheck, critical saturation points should be specified by 0 relperm - // Numerical errors and maxpermcontrast violate this even if the input has specified - // these points - if (helper.isMaster) { - if (helper.doEclipseCheck) { - for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - int minidx; - if (RelPermValues[voigtIdx][0] < RelPermValues[voigtIdx][helper.points-1]) minidx = 0; else minidx = helper.points-1; - if (RelPermValues[voigtIdx][minidx] < helper.critRelpThresh) { - RelPermValues[voigtIdx][minidx] = 0.0; - } - else { - cerr << "Minimum upscaled relperm value is " << RelPermValues[voigtIdx][minidx] << ", larger than critRelpermThresh." << endl - << "(voigtidx = " << voigtIdx << ")" << endl; - usageandexit(); - } - if (helper.upscaleBothPhases) { - if (RelPermValues2[voigtIdx][0] < RelPermValues2[voigtIdx][helper.points-1]) minidx = 0; else minidx = helper.points-1; - if (RelPermValues2[voigtIdx][minidx] < helper.critRelpThresh) { - RelPermValues2[voigtIdx][minidx] = 0.0; - } - else { - cerr << "Minimum upscaled relperm value for phase 2 is " << RelPermValues2[voigtIdx][minidx] << endl - << ", larger than critRelpermThresh.(voigtidx = " << voigtIdx << ")" << endl; - usageandexit(); - } - } - } - } - } + std::array>,2> RelPermValues; + RelPermValues[0] = helper.getRelPerm(0); + if (helper.upscaleBothPhases) + RelPermValues[1] = helper.getRelPerm(1); /********************************************************************************* * Step 9 @@ -1769,18 +1690,18 @@ try Pvalues.push_back(PvaluesVsSaturation.evaluate(SatvaluesInterp[i])); } for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - MonotCubicInterpolator RelPermVsSaturation(Satvalues, RelPermValues[voigtIdx]); - RelPermValues[voigtIdx].clear(); + MonotCubicInterpolator RelPermVsSaturation(Satvalues, RelPermValues[0][voigtIdx]); + RelPermValues[0][voigtIdx].clear(); for (int i=0; i < interpolationPoints; ++i) { - RelPermValues[voigtIdx].push_back(RelPermVsSaturation.evaluate(SatvaluesInterp[i])); + RelPermValues[0][voigtIdx].push_back(RelPermVsSaturation.evaluate(SatvaluesInterp[i])); } } if (helper.upscaleBothPhases) { for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { - MonotCubicInterpolator RelPermVsSaturation(Satvalues, RelPermValues2[voigtIdx]); - RelPermValues2[voigtIdx].clear(); + MonotCubicInterpolator RelPermVsSaturation(Satvalues, RelPermValues[1][voigtIdx]); + RelPermValues[1][voigtIdx].clear(); for (int i=0; i < interpolationPoints; ++i) { - RelPermValues2[voigtIdx].push_back(RelPermVsSaturation.evaluate(SatvaluesInterp[i])); + RelPermValues[1][voigtIdx].push_back(RelPermVsSaturation.evaluate(SatvaluesInterp[i])); } } } @@ -1798,12 +1719,12 @@ try for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) - << RelPermValues[voigtIdx][i]; + << RelPermValues[0][voigtIdx][i]; } if (helper.upscaleBothPhases) { for (int voigtIdx = 0; voigtIdx < helper.tensorElementCount; ++voigtIdx) { outputtmp << showpoint << setw(fieldwidth) << setprecision(outputprecision) - << RelPermValues2[voigtIdx][i]; + << RelPermValues[1][voigtIdx][i]; } } outputtmp << endl; @@ -1845,8 +1766,8 @@ try } for (unsigned int i=0; i < Satvalues.size(); ++i) { swofx << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i] - << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[0][i] - << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues2[0][i] + << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[0][0][i] + << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[1][0][i] << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]/100000.0 << endl; } swofx << "/" << endl; @@ -1865,8 +1786,8 @@ try } for (unsigned int i=0; i < Satvalues.size(); ++i) { swofy << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i] - << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[1][i] - << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues2[1][i] + << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[0][1][i] + << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[1][1][i] << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]/100000.0 << endl; } swofy << "/" << endl; @@ -1885,8 +1806,8 @@ try } for (unsigned int i=0; i < Satvalues.size(); ++i) { swofz << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Satvalues[i] - << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[2][i] - << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues2[2][i] + << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[0][2][i] + << showpoint << setw(fieldwidth) << setprecision(outputprecision) << RelPermValues[1][2][i] << showpoint << setw(fieldwidth) << setprecision(outputprecision) << Pvalues[i]/100000.0 << endl; } swofz << "/" << endl; diff --git a/opm/upscaling/RelPermUtils.cpp b/opm/upscaling/RelPermUtils.cpp index 2019619..d23e308 100644 --- a/opm/upscaling/RelPermUtils.cpp +++ b/opm/upscaling/RelPermUtils.cpp @@ -130,4 +130,48 @@ void RelPermUpscaleHelper::collectResults() #endif } +std::vector> RelPermUpscaleHelper::getRelPerm(int phase) const +{ + SinglePhaseUpscaler::permtensor_t zeroMatrix(3,3,(double*)0); + zero(zeroMatrix); + std::vector> RelPermValues; + if (isMaster) { + RelPermValues.resize(tensorElementCount); + // Loop over all pressure points + for (int idx=0; idx < points; ++idx) { + SinglePhaseUpscaler::permtensor_t phasePermTensor = zeroMatrix; + zero(phasePermTensor); + for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { + setVoigtValue(phasePermTensor, voigtIdx, PhasePerm[phase][idx][voigtIdx]); + } + SinglePhaseUpscaler::permtensor_t relPermTensor = zeroMatrix; + prod(phasePermTensor, permTensorInv, relPermTensor); + for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { + RelPermValues[voigtIdx].push_back(getVoigtValue(relPermTensor, voigtIdx)); + } + } + // If doEclipseCheck, critical saturation points should be specified by 0 relperm + // Numerical errors and maxpermcontrast violate this even if the input has specified + // these points + if (doEclipseCheck) { + for (int voigtIdx = 0; voigtIdx < tensorElementCount; ++voigtIdx) { + int minidx; + if (RelPermValues[voigtIdx][0] < RelPermValues[voigtIdx][points-1]) minidx = 0; else minidx = points-1; + if (RelPermValues[voigtIdx][minidx] < critRelpThresh) { + RelPermValues[voigtIdx][minidx] = 0.0; + } + else { + std::stringstream str; + str << "Minimum upscaled relperm value for phase " << phase+1 << " is " + << RelPermValues[voigtIdx][minidx] << ", larger than critRelpermThresh." << std::endl + << " (voigtidx = " << voigtIdx << ")"; + throw std::runtime_error(str.str()); + } + } + } + } + + return RelPermValues; +} + } diff --git a/opm/upscaling/RelPermUtils.hpp b/opm/upscaling/RelPermUtils.hpp index 5cbf1f2..aa48097 100644 --- a/opm/upscaling/RelPermUtils.hpp +++ b/opm/upscaling/RelPermUtils.hpp @@ -84,6 +84,12 @@ namespace Opm { //! \brief Collect results from all MPI nodes. void collectResults(); + + //! \brief Calculate relperm values from phase permeabilities. + //! \param[in] phase The phase to calculate values for (0-indexed). + //! \return The phase permeability tensor values. + //! \details First index is voigt index, second index is pressure point. + std::vector> getRelPerm(int phase) const; }; }