Suppot EQUIL + WATER + GAS

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

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())

ebos/equil/initstateequil.cc

@@ -515,6 +515,13 @@ typename PressureTable<FluidSystem,Region>::Strategy
 PressureTable<FluidSystem,Region>::
 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<MaterialLawManager, FluidSystem, Region, CellID>::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<MaterialLawManager, FluidSystem, Region, CellID>::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();
 
+        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();
 
+        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 (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 (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();
 
+        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 <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
+double PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
+materialLawCapPressGasWater() const
+{
+    return this->matLawCapPress_[this->gasPos()]
+        - this->matLawCapPress_[this->waterPos()];
+}
+
 template <class MaterialLawManager, class FluidSystem, class Region, typename CellID>
 bool PhaseSaturations<MaterialLawManager, FluidSystem, Region, CellID>::
 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 <class FluidSystem, class Region>
 void PressureTable<FluidSystem, Region>::
 equil_OWG(const Region& reg, const VSpan& span)
 {
-    // Datum depth in gas zone.  Calculate phase pressure for gas first,
+    // Datum depth in oil zone.  Calculate phase pressure for oil first,
-    // followed by oil and water if applicable.
+    // 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);
     }
 }

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.
     ///

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

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

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<TypeTag, Opm::Properties::FluidSystem>;
+    auto simulator1 = initSimulator<TypeTag>("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<TypeTag>("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;