diff --git a/opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp b/opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp index 656f2daef..3c5f86252 100644 --- a/opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp +++ b/opm/simulators/flow/equil/EquilibrationHelpers_impl.hpp @@ -446,7 +446,9 @@ operator()(double s) const fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0); fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0); fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0); - fluidState.setSaturation(phase_, s); + fluidState.setSaturation(phase_, s); + if (phase_ == FluidSystem::gasPhaseIdx) + fluidState.setSaturation(FluidSystem::waterPhaseIdx, 1.0 - s); std::array pc{0.0}; using MaterialLaw = typename MaterialLawManager::MaterialLaw; diff --git a/opm/simulators/flow/equil/InitStateEquil_impl.hpp b/opm/simulators/flow/equil/InitStateEquil_impl.hpp index 28398247c..6d1ec3958 100644 --- a/opm/simulators/flow/equil/InitStateEquil_impl.hpp +++ b/opm/simulators/flow/equil/InitStateEquil_impl.hpp @@ -578,9 +578,11 @@ deriveSaturations(const Position& x, this->setEvaluationPoint(x, reg, ptable); this->initializePhaseQuantities(); - if (ptable.waterActive()) { this->deriveWaterSat(); } if (ptable.gasActive()) { this->deriveGasSat(); } + if (ptable.waterActive()) { this->deriveWaterSat(); } + + if (this->isOverlappingTransition()) { this->fixUnphysicalTransition(); } @@ -667,36 +669,43 @@ void PhaseSaturations::deriveWa { auto& sw = this->sat_.water; - const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. - - if (this->isConstCapPress(this->waterPos())) { - // Sharp interface between phases. Can derive phase saturation - // directly from knowing where 'depth' of evaluation point is - // relative to depth of O/W contact. - sw = this->fromDepthTable(this->evalPt_.region->zwoc(), - this->waterPos(), isIncr); + const auto oilActive = this->evalPt_.ptable->oilActive(); + if (!oilActive) { + // for 2p gas+water we set the water saturation to 1.0 - sg + sw = 1.0 - this->sat_.gas; } else { - // Capillary pressure curve is non-constant, meaning there is a - // transition zone between the oil and water phases. Invert - // capillary pressure relation - // - // Pcow(Sw) = Po - Pw - // - // unless the model uses "SWATINIT". In the latter case, pick the - // saturation directly from the SWATINIT array of the pertinent - // cell. - const auto pcow = this->press_.oil - this->press_.water; + const auto isIncr = false; // dPcow/dSw <= 0 for all Sw. - if (this->swatInit_.empty()) { - sw = this->invertCapPress(pcow, this->waterPos(), isIncr); + if (this->isConstCapPress(this->waterPos())) { + // Sharp interface between phases. Can derive phase saturation + // directly from knowing where 'depth' of evaluation point is + // relative to depth of O/W contact. + sw = this->fromDepthTable(this->evalPt_.region->zwoc(), + this->waterPos(), isIncr); } else { - auto [swout, newSwatInit] = this->applySwatInit(pcow); - if (newSwatInit) + // Capillary pressure curve is non-constant, meaning there is a + // transition zone between the oil and water phases. Invert + // capillary pressure relation + // + // Pcow(Sw) = Po - Pw + // + // unless the model uses "SWATINIT". In the latter case, pick the + // saturation directly from the SWATINIT array of the pertinent + // cell. + const auto pcow = this->press_.oil - this->press_.water; + + if (this->swatInit_.empty()) { sw = this->invertCapPress(pcow, this->waterPos(), isIncr); + } else { - sw = swout; + auto [swout, newSwatInit] = this->applySwatInit(pcow); + if (newSwatInit) + sw = this->invertCapPress(pcow, this->waterPos(), isIncr); + else { + sw = swout; + } } } }