From 460f839a808f2b36cf0d6ff877bff17079a12bf5 Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Wed, 28 Sep 2016 16:19:13 +0200 Subject: [PATCH] threshold presure defaults: only consider phases which exhibit a non-zero mobility i.e., if a phase is not present in the upwind DOF, it should not be considered. this handles things analogous to the opm-simulators code. (which uses the residual saturation of the phase for the decision, but fundamentally applies the same logic.) --- applications/ebos/eclthresholdpressure.hh | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/applications/ebos/eclthresholdpressure.hh b/applications/ebos/eclthresholdpressure.hh index 30141baee..14e455e78 100644 --- a/applications/ebos/eclthresholdpressure.hh +++ b/applications/ebos/eclthresholdpressure.hh @@ -192,14 +192,22 @@ private: unsigned equilRegionOutside = elemEquilRegion_[outsideElemIdx]; if (equilRegionInside == equilRegionOutside) + // the current face is not at the boundary between EQUIL regions! continue; - // determine the maximum difference of the pressure of all phases over - // the intersection - Scalar pth = 0; + // determine the maximum difference of the pressure of any phase over the + // intersection + Scalar pth = 0.0; const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, /*timeIdx=*/0); - for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) - pth = std::max(pth, std::abs(Toolbox::value(extQuants.pressureDifference(phaseIdx)))); + for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) { + unsigned upIdx = extQuants.upstreamIndex(phaseIdx); + const auto& up = elemCtx.intensiveQuantities(upIdx, /*timeIdx=*/0); + + if (up.mobility(phaseIdx) > 0.0) { + Scalar phaseVal = Toolbox::value(extQuants.pressureDifference(phaseIdx)); + pth = std::max(pth, std::abs(phaseVal)); + } + } int offset1 = equilRegionInside*numEquilRegions_ + equilRegionOutside; int offset2 = equilRegionOutside*numEquilRegions_ + equilRegionInside;