diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 0d2d3f636..b950ac223 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -234,6 +234,8 @@ list (APPEND TEST_DATA_FILES tests/capillary_overlap.DATA tests/capillarySwatinit.DATA tests/deadfluids.DATA + tests/equil_co2store_go.DATA + tests/equil_co2store_gw.DATA tests/equil_wetgas.DATA tests/equil_liveoil.DATA tests/equil_humidwetgas.DATA diff --git a/ebos/eclequilinitializer.hh b/ebos/eclequilinitializer.hh index 79256c39b..c0901e701 100644 --- a/ebos/eclequilinitializer.hh +++ b/ebos/eclequilinitializer.hh @@ -139,12 +139,12 @@ public: if (FluidSystem::enableDissolvedGas()) fluidState.setRs(initialState.rs()[elemIdx]); - else if (Indices::gasEnabled) + else if (Indices::gasEnabled && Indices::oilEnabled) fluidState.setRs(0.0); if (FluidSystem::enableVaporizedOil()) fluidState.setRv(initialState.rv()[elemIdx]); - else if (Indices::gasEnabled) + else if (Indices::gasEnabled && Indices::oilEnabled) fluidState.setRv(0.0); if (FluidSystem::enableVaporizedWater()) diff --git a/ebos/equil/initstateequil.cc b/ebos/equil/initstateequil.cc index da93a3e44..95b4b61da 100644 --- a/ebos/equil/initstateequil.cc +++ b/ebos/equil/initstateequil.cc @@ -515,6 +515,13 @@ typename PressureTable::Strategy PressureTable:: selectEquilibrationStrategy(const Region& reg) const { + if (!this->oilActive()) { + if (reg.datum() > reg.zwoc()) { // Datum in water zone + return &PressureTable::equil_WOG; + } + return &PressureTable::equil_GOW; + } + if (reg.datum() > reg.zwoc()) { // Datum in water zone return &PressureTable::equil_WOG; } @@ -636,12 +643,14 @@ void PhaseSaturations::deriveGa auto& sg = this->sat_.gas; const auto isIncr = true; // dPcgo/dSg >= 0 for all Sg. + const auto oilActive = this->evalPt_.ptable->oilActive(); if (this->isConstCapPress(this->gasPos())) { // Sharp interface between phases. Can derive phase saturation // directly from knowing where 'depth' of evaluation point is // relative to depth of O/G contact. - sg = this->fromDepthTable(this->evalPt_.region->zgoc(), + const auto gas_contact = oilActive? this->evalPt_.region->zgoc() : this->evalPt_.region->zwoc(); + sg = this->fromDepthTable(gas_contact, this->gasPos(), isIncr); } else { @@ -652,8 +661,8 @@ void PhaseSaturations::deriveGa // Pcgo(Sg) = Pg - Po // // Note that Pcgo is defined to be (Pg - Po), not (Po - Pg). - const auto pcgo = this->press_.gas - this->press_.oil; - + const auto pw = oilActive? this->press_.oil : this->press_.water; + const auto pcgo = this->press_.gas - pw; sg = this->invertCapPress(pcgo, this->gasPos(), isIncr); } } @@ -732,64 +741,62 @@ accountForScaledSaturations() { const auto gasActive = this->evalPt_.ptable->gasActive(); const auto watActive = this->evalPt_.ptable->waterActive(); + const auto oilActive = this->evalPt_.ptable->oilActive(); + + auto sg = gasActive? this->sat_.gas : 0.0; + auto sw = watActive? this->sat_.water : 0.0; + auto so = oilActive? this->sat_.oil : 0.0; + + this->fluidState_.setSaturation(this->waterPos(), sw); + this->fluidState_.setSaturation(this->oilPos(), so); + this->fluidState_.setSaturation(this->gasPos(), sg); const auto& scaledDrainageInfo = this->matLawMgr_ .oilWaterScaledEpsInfoDrainage(this->evalPt_.position->cell); - const auto sg = this->sat_.gas; - const auto sw = this->sat_.water; - - { - auto so = 1.0; - - if (watActive) { - const auto swu = scaledDrainageInfo.Swu; - so -= swu; - - this->fluidState_.setSaturation(this->waterPos(), swu); - } - - if (gasActive) { - const auto sgu = scaledDrainageInfo.Sgu; - so -= sgu; - - this->fluidState_.setSaturation(this->gasPos(), sgu); - } - - this->fluidState_.setSaturation(this->oilPos(), so); - } - const auto thresholdSat = 1.0e-6; if (watActive && ((sw + thresholdSat) > scaledDrainageInfo.Swu)) { // Water saturation exceeds maximum possible value. Reset oil phase // pressure to that which corresponds to maximum possible water // saturation value. this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swu); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swu); + } else if (gasActive) { + this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swu); + } + sw = scaledDrainageInfo.Swu; this->computeMaterialLawCapPress(); - // Pcow = Po - Pw => Po = Pw + Pcow - this->press_.oil = this->press_.water + this->materialLawCapPressOilWater(); + if (oilActive) { + // Pcow = Po - Pw => Po = Pw + Pcow + this->press_.oil = this->press_.water + this->materialLawCapPressOilWater(); + } else { + // Pcgw = Pg - Pw => Pg = Pw + Pcgw + this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); + } + } - else if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) { + if (gasActive && ((sg + thresholdSat) > scaledDrainageInfo.Sgu)) { // Gas saturation exceeds maximum possible value. Reset oil phase // pressure to that which corresponds to maximum possible gas // saturation value. this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgu); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgu); + } else if (watActive) { + this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgu); + } + sg = scaledDrainageInfo.Sgu; this->computeMaterialLawCapPress(); - // Pcgo = Pg - Po => Po = Pg - Pcgo - this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); - } - - if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) { - // Gas saturation less than minimum possible value in cell. Reset - // gas phase pressure to that which corresponds to minimum possible - // gas saturation. - this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl); - this->computeMaterialLawCapPress(); - - // Pcgo = Pg - Po => Pg = Po + Pcgo - this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil(); + if (oilActive) { + // Pcgo = Pg - Po => Po = Pg - Pcgo + this->press_.oil = this->press_.gas - this->materialLawCapPressGasOil(); + } else { + // Pcgw = Pg - Pw => Pw = Pg - Pcgw + this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); + } } if (watActive && ((sw - thresholdSat) < scaledDrainageInfo.Swl)) { @@ -797,10 +804,43 @@ accountForScaledSaturations() // water phase pressure to that which corresponds to minimum // possible water saturation value. this->fluidState_.setSaturation(this->waterPos(), scaledDrainageInfo.Swl); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sw - scaledDrainageInfo.Swl); + } else if (gasActive) { + this->fluidState_.setSaturation(this->gasPos(), sg + sw - scaledDrainageInfo.Swl); + } + sw = scaledDrainageInfo.Swl; this->computeMaterialLawCapPress(); - // Pcwo = Po - Pw => Pw = Po - Pcow - this->press_.water = this->press_.oil - this->materialLawCapPressOilWater(); + if (oilActive) { + // Pcwo = Po - Pw => Pw = Po - Pcow + this->press_.water = this->press_.oil - this->materialLawCapPressOilWater(); + } else { + // Pcgw = Pg - Pw => Pw = Pg - Pcgw + this->press_.water = this->press_.gas - this->materialLawCapPressGasWater(); + } + } + + if (gasActive && ((sg - thresholdSat) < scaledDrainageInfo.Sgl)) { + // Gas saturation less than minimum possible value in cell. Reset + // gas phase pressure to that which corresponds to minimum possible + // gas saturation. + this->fluidState_.setSaturation(this->gasPos(), scaledDrainageInfo.Sgl); + if (oilActive) { + this->fluidState_.setSaturation(this->oilPos(), so + sg - scaledDrainageInfo.Sgl); + } else if (watActive) { + this->fluidState_.setSaturation(this->waterPos(), sw + sg - scaledDrainageInfo.Sgl); + } + sg = scaledDrainageInfo.Sgl; + this->computeMaterialLawCapPress(); + + if (oilActive) { + // Pcgo = Pg - Po => Pg = Po + Pcgo + this->press_.gas = this->press_.oil + this->materialLawCapPressGasOil(); + } else { + // Pcgw = Pg - Pw => Pg = Pw + Pcgw + this->press_.gas = this->press_.water + this->materialLawCapPressGasWater(); + } } } @@ -847,6 +887,14 @@ materialLawCapPressOilWater() const - this->matLawCapPress_[this->waterPos()]; } +template +double PhaseSaturations:: +materialLawCapPressGasWater() const +{ + return this->matLawCapPress_[this->gasPos()] + - this->matLawCapPress_[this->waterPos()]; +} + template bool PhaseSaturations:: isConstCapPress(const PhaseIdx phaseIdx) const @@ -1081,7 +1129,6 @@ equil_GOW(const Region& reg, const VSpan& span) reg.zgoc(), this->gas(reg.zgoc()) - reg.pcgoGoc() }; - this->makeOilPressure(ic, reg, span); } @@ -1107,8 +1154,8 @@ template void PressureTable:: equil_OWG(const Region& reg, const VSpan& span) { - // Datum depth in gas zone. Calculate phase pressure for gas first, - // followed by oil and water if applicable. + // Datum depth in oil zone. Calculate phase pressure for oil first, + // followed by gas and water if applicable. if (! this->oilActive()) { throw std::invalid_argument { @@ -1141,7 +1188,6 @@ equil_OWG(const Region& reg, const VSpan& span) reg.zgoc(), this->oil(reg.zgoc()) + reg.pcgoGoc() }; - this->makeGasPressure(ic, reg, span); } } diff --git a/ebos/equil/initstateequil.hh b/ebos/equil/initstateequil.hh index af995160c..8b7096a09 100644 --- a/ebos/equil/initstateequil.hh +++ b/ebos/equil/initstateequil.hh @@ -574,6 +574,10 @@ private: /// fluid state. double materialLawCapPressOilWater() const; + /// Extract gas/water capillary pressure value (Pg - Pw) from current + /// fluid state. + double materialLawCapPressGasWater() const; + /// Predicate for whether specific phase has constant capillary pressure /// curve in current cell. /// diff --git a/tests/equil_co2store_go.DATA b/tests/equil_co2store_go.DATA new file mode 100644 index 000000000..2ff8efcd0 --- /dev/null +++ b/tests/equil_co2store_go.DATA @@ -0,0 +1,106 @@ +NOECHO + +RUNSPEC ====== + +OIL +GAS +CO2STORE + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 20 +/ + +WELLDIMS + 30 10 2 30 / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +20*5.0 +/ + + +PORO +20*0.2 +/ + + +PERMZ + 20*1.0 +/ + +PERMY +20*100.0 +/ + +PERMX +20*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + +SGOF +0 0 1 0.1 +0.9 1 0 0.5 +/ + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +SOLUTION ====== + +EQUIL +100 150 1000 1* 50 0.2 1* 1* 0 +/ + +RPTSOL +'PRES' 'PGAS' 'SOIL' 'SGAS' 'RS' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +TSTEP +1 / + +RPTSCHED +'PRES' 'PGAS' 'SOIL' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/equil_co2store_gw.DATA b/tests/equil_co2store_gw.DATA new file mode 100644 index 000000000..6122549ed --- /dev/null +++ b/tests/equil_co2store_gw.DATA @@ -0,0 +1,111 @@ +NOECHO + +RUNSPEC ====== + +WATER +GAS +CO2STORE + +TABDIMS + 1 1 40 20 1 20 / + +DIMENS +1 1 20 +/ + +WELLDIMS + 30 10 2 30 / + +START + 1 'JAN' 1990 / + +NSTACK + 25 / + +EQLDIMS +-- NTEQUL + 1 / + + +FMTOUT +FMTIN + +GRID ====== + +DXV +1.0 +/ + +DYV +1.0 +/ + +DZV +20*5.0 +/ + + +PORO +20*0.2 +/ + + +PERMZ + 20*1.0 +/ + +PERMY +20*100.0 +/ + +PERMX +20*100.0 +/ + +BOX + 1 1 1 1 1 1 / + +TOPS +0.0 +/ + +PROPS ====== + +SWFN +0.1 0 0.5 +1 1 0.1 +/ + +SGFN +0 0 0.0 +0.9 1 0.0 +/ + +ROCK +--RefPres Comp + 1. 5.0E-5 / + +SOLUTION ====== + +EQUIL +100 150 50 0.2 1000 1* 1* 1* 0 +/ + +RPTSOL +'PRES' 'PGAS' 'PWAT' 'SWAT' 'SGAS' 'RSW' 'RESTART=2' / + +SUMMARY ====== +RUNSUM + +SEPARATE + +SCHEDULE ====== + +TSTEP +1 / + +RPTSCHED +'PRES' 'PGAS' 'PWAT' 'SWAT' 'SGAS' 'RSW' 'RESTART=3' 'NEWTON=2' / + + +END diff --git a/tests/test_equil.cc b/tests/test_equil.cc index a29196718..b0c99612c 100644 --- a/tests/test_equil.cc +++ b/tests/test_equil.cc @@ -790,6 +790,49 @@ BOOST_AUTO_TEST_CASE(DeckWithLiveOil) } } +BOOST_AUTO_TEST_CASE(DeckWithCO2STORE) +{ + using TypeTag = Opm::Properties::TTag::TestEquilTypeTag; + using FluidSystem = Opm::GetPropType; + auto simulator1 = initSimulator("equil_co2store_go.DATA"); + EquilFixture::Initializer comp_go(*simulator1->problem().materialLawManager(), + simulator1->vanguard().eclState(), + simulator1->vanguard().grid(), + simulator1->vanguard().gridView(), + simulator1->vanguard().cartesianMapper(), 9.80665); + + auto simulator2 = initSimulator("equil_co2store_gw.DATA"); + EquilFixture::Initializer comp_gw(*simulator2->problem().materialLawManager(), + simulator2->vanguard().eclState(), + simulator2->vanguard().grid(), + simulator2->vanguard().gridView(), + simulator2->vanguard().cartesianMapper(), 9.80665); + + Opm::GridManager gm(simulator2->vanguard().eclState().getInputGrid()); + const UnstructuredGrid& grid = *(gm.c_grid()); + + const double reltol = 1.0e-5; + const auto& pressures_go = comp_go.press(); + BOOST_REQUIRE_EQUAL(pressures_go.size(), 3U); + BOOST_REQUIRE_EQUAL(int(pressures_go[0].size()), grid.number_of_cells); + + const auto& pressures_gw = comp_gw.press(); + BOOST_REQUIRE_EQUAL(pressures_gw.size(), 3U); + BOOST_REQUIRE_EQUAL(int(pressures_gw[0].size()), grid.number_of_cells); + + const auto& sats_go = comp_go.saturation(); + const auto& sats_gw = comp_gw.saturation(); + + for (int i = 0; i < grid.number_of_cells; ++i) { + BOOST_CHECK_CLOSE(pressures_go[FluidSystem::gasPhaseIdx][i], pressures_gw[FluidSystem::gasPhaseIdx][i], reltol); + BOOST_CHECK_CLOSE(pressures_go[FluidSystem::oilPhaseIdx][i], pressures_gw[FluidSystem::waterPhaseIdx][i], reltol); + + BOOST_CHECK_CLOSE(sats_go[FluidSystem::gasPhaseIdx][i], sats_gw[FluidSystem::gasPhaseIdx][i], reltol); + BOOST_CHECK_CLOSE(sats_go[FluidSystem::oilPhaseIdx][i], sats_gw[FluidSystem::waterPhaseIdx][i], reltol); + } +} + + BOOST_AUTO_TEST_CASE(DeckWithWetGas) { using TypeTag = Opm::Properties::TTag::TestEquilTypeTag;