adding primary variable combination/switching relevant for combined vaporization of water and oil into gas phase

opm/models/blackoil/blackoilbrinemodules.hh

@@ -467,7 +467,7 @@ public:
                                   unsigned timeIdx)
     {
         if (enableSaltPrecipitation) {
-            const Evaluation porosityFactor  = 1.0 - saltSaturation(); //phi/phi_0
+            const Evaluation porosityFactor  = min(1.0 - saltSaturation(), 1.0); //phi/phi_0
 
             const auto& permfactTable = BrineModule::permfactTable(elemCtx, dofIdx, timeIdx);
 

opm/models/blackoil/blackoilintensivequantities.hh

@@ -149,7 +149,7 @@ public:
         if (waterEnabled) {
             if (priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p) {
                 Sw = 1.0;
-            } else if (priVars.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg) {
+            } else if (priVars.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg && priVars.primaryVarsMeaning() != PrimaryVariables::Rvw_pg_Rv) {
                 Sw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
             }
         }
@@ -167,8 +167,14 @@ public:
                 // deal with solvent
                 if (enableSolvent)
                     Sg -= priVars.makeEvaluation(Indices::solventSaturationIdx, timeIdx);
-            }
+
-            else
+            } else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_po_Sg) {
+                 // -> oil-gas case
+                 Sg = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            } else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
+                 // -> gas case
+                 Sg = 1.0;
+            } else
             {
                 assert(priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs);
                 // -> oil-water case
@@ -205,7 +211,7 @@ public:
         MaterialLaw::capillaryPressures(pC, materialParams, fluidState_);
 
         //oil is the reference phase for pressure 
-        if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
+        if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv || priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
             const Evaluation& pg = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
             for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx)
                 if (FluidSystem::phaseIsActive(phaseIdx))
@@ -276,6 +282,42 @@ public:
             // The switching variable is the water-gas ratio Rvw
             const auto& Rvw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
             fluidState_.setRvw(Rvw);
+
+            if (FluidSystem::enableVaporizedOil()) {
+                Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
+                const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+                fluidState_.setRv(min(RvMax, RvSat));
+            }
+            else if (compositionSwitchEnabled)
+                fluidState_.setRv(0.0);
+        }
+        else if (priVars.primaryVarsMeaning() == PrimaryVariables::Rvw_pg_Rv) {
+            // The switching variable is the water-gas ratio Rvw
+            const auto& Rvw = priVars.makeEvaluation(Indices::waterSaturationIdx, timeIdx);
+            fluidState_.setRvw(Rvw);
+
+            const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
+            fluidState_.setRv(Rv);
+
+            if (FluidSystem::enableDissolvedGas()) {
+                // the oil phase is not present, but we need to compute its "composition" for
+                // the gravity correction anyway
+                Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
+                const auto& RsSat = enableExtbo ? asImp_().rs() :
+                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                                                            oilPhaseIdx,
+                                                            pvtRegionIdx,
+                                                            SoMax);
+
+                fluidState_.setRs(min(RsMax, RsSat));
+            } 
+            else {
+                fluidState_.setRs(0.0);
+            }
         }
         else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_po_Rs) {
             // if the switching variable is the mole fraction of the gas component in the
@@ -299,6 +341,13 @@ public:
             }
             else
                 fluidState_.setRv(0.0);
+
+            if (FluidSystem::enableVaporizedWater()) {
+                const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx);
+                fluidState_.setRvw(RvwSat);
+            }
         }
         else if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
             const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
@@ -318,6 +367,13 @@ public:
             } else {
                 fluidState_.setRs(0.0);
             }
+
+            if (FluidSystem::enableVaporizedWater()) {
+                const Evaluation& RvwSat = FluidSystem::saturatedVaporizationFactor(fluidState_,
+                                                            gasPhaseIdx,
+                                                            pvtRegionIdx);
+                fluidState_.setRvw(RvwSat);
+            }
         } else {
             assert(priVars.primaryVarsMeaning() == PrimaryVariables::OnePhase_p);
         }
@@ -364,6 +420,12 @@ public:
                     fluidState_.Rv() *
                     FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
             }
+            if (FluidSystem::enableVaporizedWater()) {
+                rho +=
+                    fluidState_.invB(gasPhaseIdx) *
+                    fluidState_.Rvw() *
+                    FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
+            }
             fluidState_.setDensity(gasPhaseIdx, rho);
         }
 
@@ -409,6 +471,12 @@ public:
           porosity_ += - biofilm_ - calcite_;
         }
 
+        // deal with salt-precipitation
+        if (enableSaltPrecipitation && priVars.primaryVarsMeaningBrine() == PrimaryVariables::Sp) {
+            Evaluation Sp = priVars.makeEvaluation(Indices::saltConcentrationIdx, timeIdx);
+            porosity_ *= (1.0 - Sp);
+        }
+
         asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().zPvtUpdate_();
         asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -264,12 +264,15 @@ protected:
         Scalar deltaSg = 0.0;
         Scalar deltaSs = 0.0;
 
-        if (Indices::waterEnabled && FluidSystem::numActivePhases() > 1 && currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg) {
+        if (Indices::waterEnabled && FluidSystem::numActivePhases() > 1
+            && currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg 
+            && currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_pg_Rv) {
             deltaSw = update[Indices::waterSaturationIdx];
             deltaSo = -deltaSw;
         }
 
-        if (compositionSwitchEnabled && currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg) {
+        if (compositionSwitchEnabled && (currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg
+            || currentValue.primaryVarsMeaning() == PrimaryVariables::Rvw_po_Sg)) {
             deltaSg = update[Indices::compositionSwitchIdx];
             deltaSo -= deltaSg;
         }
@@ -306,7 +309,8 @@ protected:
             }
             // water saturation delta
             else if (pvIdx == Indices::waterSaturationIdx)
-                if (currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg)
+                if (currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_po_Sg 
+                    && currentValue.primaryVarsMeaning() != PrimaryVariables::Rvw_pg_Rv)
                     delta *= satAlpha;
                 else{
                     //Ensure Rvw factor does not become negative
@@ -319,7 +323,8 @@ protected:
                 // interpretation since it can represent Sg, Rs or Rv. For now, we only
                 // limit saturation deltas and ensure that the R factors do not become
                 // negative.
-                if (currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg)
+                if (currentValue.primaryVarsMeaning() == PrimaryVariables::Sw_po_Sg 
+                    || currentValue.primaryVarsMeaning() == PrimaryVariables::Rvw_po_Sg)
                     delta *= satAlpha;
                 else {
                     if (delta > currentValue[Indices::compositionSwitchIdx])

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -129,7 +129,7 @@ public:
         Sw_po_Rs, // water + oil case
         Sw_pg_Rv, // water + gas case
         Rvw_po_Sg, // gas + oil  case
-
+        Rvw_pg_Rv, //gas only
         OnePhase_p, // onephase case
     };
 
@@ -308,7 +308,10 @@ public:
         else if ((gasPresent && oilPresent) || (onlyWater && FluidSystem::phaseIsActive(oilPhaseIdx))) {
             // gas and oil: both hydrocarbon phases are in equilibrium (i.e., saturated
             // with the "protagonist" component of the other phase.)
-            primaryVarsMeaning_ = Sw_po_Sg;
+            if (FluidSystem::enableVaporizedWater() && !waterPresent)
+                primaryVarsMeaning_ = Rvw_po_Sg; 
+            else
+                primaryVarsMeaning_ = Sw_po_Sg;
         }
         else if (oilPresent) {
             // only oil: if dissolved gas is enabled, we need to consider the oil phase
@@ -380,6 +383,21 @@ public:
             }
             if( compositionSwitchEnabled )
                 (*this)[compositionSwitchIdx] = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
+        }
+         else if (primaryVarsMeaning() == Rvw_pg_Rv && FluidSystem::enableVaporizedWater()) {
+            const auto& Rvw = BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
+            const auto& Rv = BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fluidState, pvtRegionIdx_);
+            if (waterEnabled)
+                (*this)[waterSaturationIdx] = Rvw; //waterSaturationIdx becomes a switching idx
+            if (gasEnabled && waterEnabled && !oilEnabled) {
+                //-> water-gas system
+                (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(gasPhaseIdx));
+            }
+            else if (oilEnabled) {
+                (*this)[pressureSwitchIdx] = FsToolbox::value(fluidState.pressure(oilPhaseIdx));
+            }
+            if( compositionSwitchEnabled )
+                (*this)[compositionSwitchIdx] = Rv;
         }
         else {
             assert(primaryVarsMeaning() == Sw_pg_Rv);
@@ -432,7 +450,7 @@ public:
             return false;
         }
         Scalar Sw = 0.0;
-        if (waterEnabled && primaryVarsMeaning() != Rvw_po_Sg)
+        if (waterEnabled && primaryVarsMeaning() != Rvw_po_Sg && primaryVarsMeaning() != Rvw_pg_Rv )
             Sw = (*this)[Indices::waterSaturationIdx];
 
         if (enableSaltPrecipitation) {
@@ -474,14 +492,14 @@ public:
                 Sg = (*this)[Indices::compositionSwitchIdx];
 
             Scalar So = 1.0 - Sw - Sg - solventSaturation_();
-
+            Scalar So3 = 1.0 - Sg - solventSaturation_();
             //water disappears 
-            if( Sw < -eps && FluidSystem::enableVaporizedWater()) {
+            if(Sw < -eps && So3 > 0.0 && Sg > 0.0  && FluidSystem::enableVaporizedWater()) {
                 Scalar po = (*this)[Indices::pressureSwitchIdx]; 
                 Scalar T = asImp_().temperature_();
                 Scalar pC[numPhases] = { 0.0 };
                 const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
-                computeCapillaryPressures_(pC, So, Sg + solventSaturation_(), Sw, matParams);
+                computeCapillaryPressures_(pC, So3, Sg + solventSaturation_(), /*Sw=*/ 0.0, matParams);
                 Scalar pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
                 Scalar RvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
                                                                                       T,
@@ -491,6 +509,36 @@ public:
 
                 return true;
             }
+            //water and oil disappears 
+            if(Sw < -eps && So3 <-eps && Sg > 0.0 && FluidSystem::enableVaporizedWater() && FluidSystem::enableVaporizedOil()) {
+                Scalar po = (*this)[Indices::pressureSwitchIdx]; 
+                Scalar T = asImp_().temperature_();
+                Scalar pC[numPhases] = { 0.0 };
+                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                computeCapillaryPressures_(pC, /*So*/ 0.0, Sg + solventSaturation_(), /*Sw=*/ 0.0, matParams);
+                Scalar pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
+                Scalar RvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
+                                                                                      T,
+                                                                                      pg);
+
+                Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
+                Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
+                Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
+                                                             pg,
+                                                             zFraction_())
+                             : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                                                    T,
+                                                                                    pg,
+                                                                                    Scalar(0),
+                                                                                    SoMax); 
+                setPrimaryVarsMeaning(Rvw_pg_Rv);
+                (*this)[Indices::pressureSwitchIdx] = pg;
+                (*this)[Indices::waterSaturationIdx] = RvwSat; //primary variable becomes Rvw 
+                if (compositionSwitchEnabled)
+                    (*this)[Indices::compositionSwitchIdx] = std::min(RvMax, RvSat);
+
+                return true;
+            }
 
             Scalar So2 = 1.0 - Sw - solventSaturation_();
             if (Sg < -eps && So2 > 0.0 && FluidSystem::enableDissolvedGas()) {
@@ -620,51 +668,107 @@ public:
                                                                                    T,
                                                                                    pg); 
             Scalar Rvw = (*this)[Indices::waterSaturationIdx];
-            if (Rvw >  RvwSat*(1.0 + eps)) {
+            if (Rvw > RvwSat*(1.0 + eps)) {
                 // water phase appears
                 setPrimaryVarsMeaning(Sw_po_Sg);
                 (*this)[Indices::waterSaturationIdx] = 0.0; // water saturation
 
                 return true;
             }
+            if (So < -eps && FluidSystem::enableVaporizedOil()) {
+                //oil phase dissappears
+                computeCapillaryPressures_(pC, /*So=*/ 0.0, Sg + solventSaturation_(), /*Sw=*/ 0.0, matParams);
+                pg = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
+                 Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
+                Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
+                Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
+                                                             pg,
+                                                             zFraction_())
+                             : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                                                    T,
+                                                                                    pg,
+                                                                                    Scalar(0),
+                                                                                    SoMax);
+                setPrimaryVarsMeaning(Rvw_pg_Rv);
+                (*this)[Indices::pressureSwitchIdx] = pg;
+                if (compositionSwitchEnabled)
+                    (*this)[Indices::compositionSwitchIdx] = std::min(RvMax, RvSat);
+
+                return true;
+            }
+            return false;
+        }
+        else if (primaryVarsMeaning() == Rvw_pg_Rv) {
+            //only gas phase
+            Scalar pg = (*this)[Indices::pressureSwitchIdx]; 
+            Scalar T = asImp_().temperature_();
+            Scalar RvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
+                                                                                   T,
+                                                                                   pg); 
+            Scalar Rvw = (*this)[Indices::waterSaturationIdx];
+            if (Rvw > RvwSat*(1.0 + eps)) {
+                // water phase appears
+                setPrimaryVarsMeaning(Sw_pg_Rv);
+                (*this)[Indices::waterSaturationIdx] = 0.0; // water saturation
+
+                return true;
+            }
+
+            Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
+            Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
+            Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
+                                                         pg,
+                                                         zFraction_())
+                         : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                                                T,
+                                                                                pg,
+                                                                                /*So=*/Scalar(0.0),
+                                                                                SoMax);
+
+            Scalar Rv = (*this)[Indices::compositionSwitchIdx];
+            if (Rv > std::min(RvMax, RvSat*(1.0 + eps))) {
+                // switch to phase equilibrium mode because the oil phase appears. here
+                // we also need the capillary pressures to calculate the oil phase
+                // pressure using the gas phase pressure
+                Scalar pC[numPhases] = { 0.0 };
+                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                computeCapillaryPressures_(pC,
+                                           /*So=*/0.0,
+                                           /*Sg=*/1.0,
+                                           /*Sw=*/0.0,
+                                           matParams);
+                Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
+
+                setPrimaryVarsMeaning(Rvw_po_Sg);
+                (*this)[Indices::pressureSwitchIdx] = po;
+                (*this)[Indices::compositionSwitchIdx] = 1.0; // hydrocarbon gas saturation
+
+                return true;
+            }
             return false;
         }    
         else {
             assert(primaryVarsMeaning() == Sw_pg_Rv);
             assert(compositionSwitchEnabled);
 
+            Scalar T = asImp_().temperature_();
             Scalar pg = (*this)[Indices::pressureSwitchIdx];
             Scalar Sg = 1.0 - Sw - solventSaturation_();
 
-            // special case for cells with almost only water
+            //water disappears 
-            if (Sw >= thresholdWaterFilledCell) {
+            if(Sw < -eps  && FluidSystem::enableVaporizedWater()) {
-                // switch to phase equilibrium mode because the hydrocarbon gas phase
+                Scalar RvwSat = FluidSystem::gasPvt().saturatedWaterVaporizationFactor(pvtRegionIdx_,
-                // disappears. here we need the capillary pressures to calculate the oil
+                                                                                      T,
-                // phase pressure using the gas phase pressure
+                                                                                      pg);
-                Scalar pC[numPhases] = { 0.0 };
+                setPrimaryVarsMeaning(Rvw_pg_Rv);
-                const MaterialLawParams& matParams = problem.materialLawParams(globalDofIdx);
+                (*this)[Indices::waterSaturationIdx] = RvwSat; //primary variable becomes Rvw 
-                computeCapillaryPressures_(pC,
-                                           /*So=*/0.0,
-                                           /*Sg=*/Sg + solventSaturation_(),
-                                           Sw,
-                                           matParams);
-                Scalar po = pg + (pC[oilPhaseIdx] - pC[gasPhaseIdx]);
-
-                setPrimaryVarsMeaning(Sw_po_Sg);
-                if (waterEnabled)
-                    (*this)[Indices::waterSaturationIdx] = 1.0;
-
-                (*this)[Indices::pressureSwitchIdx] = po;
-                (*this)[Indices::compositionSwitchIdx] = 0.0; // hydrocarbon gas saturation
 
                 return true;
             }
-
             // Only the gas and the water phases are present. The oil phase appears as
             // soon as more of the oil component is present in the hydrocarbon gas phase
             // than what saturated gas contains. Note that we use the blackoil specific
             // low-level PVT objects here for performance reasons.
-            Scalar T = asImp_().temperature_();
             Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
             Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
             Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),