From 83a9f6ffce1f93eecee71fc753031fbd68c8af33 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 28 Sep 2016 16:38:41 +0200 Subject: [PATCH] threshold pressure defaults: calculate the pressure difference like when computing the fluxes this should not change the value of the result at all (because the total delta which is added to the phase pressures stays identical), but it should be less confusing when comparing this with the code that calculates the gravity correction term in the flux calculation. --- opm/simulators/thresholdPressures.hpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/opm/simulators/thresholdPressures.hpp b/opm/simulators/thresholdPressures.hpp index 9fdcb1ddc..6c7b5ff09 100644 --- a/opm/simulators/thresholdPressures.hpp +++ b/opm/simulators/thresholdPressures.hpp @@ -278,7 +278,6 @@ void computeMaxDp(std::map, double>& maxDp, for (int phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { const double z1 = UgGridHelpers::cellCenterDepth(grid, c1); const double z2 = UgGridHelpers::cellCenterDepth(grid, c2); - const double zAvg = (z1 + z2)/2; // average depth const double rhoAvg = (rho[phaseIdx][c1] + rho[phaseIdx][c2])/2; @@ -289,8 +288,8 @@ void computeMaxDp(std::map, double>& maxDp, const double sResid2 = minSat[numPhases*c2 + phaseIdx]; // compute gravity corrected pressure potentials at the average depth - const double p1 = phasePressure[phaseIdx][c1] + rhoAvg*gravity*(zAvg - z1); - const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(zAvg - z2); + const double p1 = phasePressure[phaseIdx][c1]; + const double p2 = phasePressure[phaseIdx][c2] + rhoAvg*gravity*(z1 - z2); if ((p1 > p2 && s1 > sResid1) || (p2 > p1 && s2 > sResid2)) maxDp[barrierId] = std::max(maxDp[barrierId], std::abs(p1 - p2));