extending drsdtcon with regimes, and option for GAS/WATER

opm/models/blackoil/blackoilconvectivemixingmodule.hh

@@ -0,0 +1,269 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 2 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+
+  Consult the COPYING file in the top-level source directory of this
+  module for the precise wording of the license and the list of
+  copyright holders.
+*/
+/*!
+ * \file
+ *
+ * \brief Classes required for molecular diffusion.
+ */
+#ifndef EWOMS_CONVECTIVEMIXING_MODULE_HH
+#define EWOMS_CONVECTIVEMIXING_MODULE_HH
+
+#include <opm/models/discretization/common/fvbaseproperties.hh>
+#include <opm/input/eclipse/Schedule/OilVaporizationProperties.hpp>
+#include <opm/input/eclipse/Schedule/Schedule.hpp>
+#include <opm/material/common/Valgrind.hpp>
+
+#include <dune/common/fvector.hh>
+
+#include <stdexcept>
+#include <iostream>
+
+namespace Opm {
+
+/*!
+ * \copydoc Opm::BlackOilConvectiveMixingModule
+ * \brief Provides the requiered methods for dynamic convective mixing.
+ */
+template <class TypeTag>
+class BlackOilConvectiveMixingModule
+{
+    using Scalar = GetPropType<TypeTag, Properties::Scalar>;
+    using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>;
+    using RateVector = GetPropType<TypeTag, Properties::RateVector>;
+    using Evaluation = GetPropType<TypeTag, Properties::Evaluation>;
+    using Indices = GetPropType<TypeTag, Properties::Indices>;
+    using IntensiveQuantities = GetPropType<TypeTag, Properties::IntensiveQuantities>;
+    using GridView = GetPropType<TypeTag, Properties::GridView>;
+
+    enum { conti0EqIdx = Indices::conti0EqIdx };
+    enum { dimWorld = GridView::dimensionworld };
+
+public:
+
+    struct ConvectiveMixingModuleParam
+    {
+        bool active_;
+        std::vector<Scalar> Xhi_;
+        std::vector<Scalar> Psi_;
+    };
+
+    #if HAVE_ECL_INPUT
+    static void beginEpisode(const EclipseState& eclState, const Schedule& schedule, const int episodeIdx, ConvectiveMixingModuleParam& info)
+    {
+        // check that Xhi and Psi didn't change
+        std::size_t numRegions = eclState.runspec().tabdims().getNumPVTTables();
+        const auto& control = schedule[episodeIdx].oilvap();
+        if (!info.active_) {
+            return;
+        }
+        if (info.Xhi_.empty()) {
+            info.Xhi_.resize(numRegions); 
+            info.Psi_.resize(numRegions); 
+        }
+        for (size_t i = 0; i < numRegions; ++i ) {
+            info.Xhi_[i] = control.getMaxDRSDT(i);
+            info.Psi_[i] = control.getPsi(i);
+        }
+    }
+    #endif
+
+    template <class Context>
+    static void addConvectiveMixingFlux(RateVector& flux, 
+                                        const Context& elemCtx,
+                                        unsigned scvfIdx, 
+                                        unsigned timeIdx)
+    {
+
+        // need for dary flux calculation
+        const auto& problem = elemCtx.problem();
+        const auto& stencil = elemCtx.stencil(timeIdx);
+        const auto& scvf = stencil.interiorFace(scvfIdx);
+
+        unsigned interiorDofIdx = scvf.interiorIndex();
+        unsigned exteriorDofIdx = scvf.exteriorIndex();
+        assert(interiorDofIdx != exteriorDofIdx);
+
+        const auto& globalIndexIn = stencil.globalSpaceIndex(interiorDofIdx);
+        const auto& globalIndexEx = stencil.globalSpaceIndex(exteriorDofIdx);
+        Scalar trans = problem.transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
+        Scalar faceArea = scvf.area();
+        const Scalar g = problem.gravity()[dimWorld - 1];
+        const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
+        const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
+        const Scalar zIn = problem.dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
+        const Scalar zEx = problem.dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
+        const Scalar distZ = zIn - zEx;
+        addConvectiveMixingFlux(flux,
+                                intQuantsIn,
+                                intQuantsEx,
+                                globalIndexIn,
+                                globalIndexEx,
+                                distZ * g,
+                                trans,
+                                faceArea,
+                                problem.moduleParams().convectiveMixingModuleParam);
+    }
+
+
+
+
+
+    /*!
+     * \brief Adds the diffusive mass flux flux to the flux vector over a flux
+     *        integration point.
+      */
+    static void addConvectiveMixingFlux(RateVector& flux,
+                            const IntensiveQuantities& intQuantsIn,
+                            const IntensiveQuantities& intQuantsEx,
+                            const unsigned globalIndexIn,
+                            const unsigned globalIndexEx,
+                            const Scalar distZg, 
+                            const Scalar trans,
+                            const Scalar faceArea,
+                            const ConvectiveMixingModuleParam& info)
+    {
+
+
+        if (!info.active_) {
+            return;
+        }
+
+        const auto& liquidPhaseIdx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
+            FluidSystem::waterPhaseIdx : 
+            FluidSystem::oilPhaseIdx;
+        const Evaluation SoMax = 0.0;
+
+        //interiour
+        const Scalar rs_zero_in = 0.0;
+        const auto& t_in = Opm::getValue(intQuantsIn.fluidState().temperature(liquidPhaseIdx));
+        const auto& p_in = Opm::getValue(intQuantsIn.fluidState().pressure(liquidPhaseIdx));
+        const auto& rssat_in = FluidSystem::saturatedDissolutionFactor(intQuantsIn.fluidState(),
+                                                            liquidPhaseIdx,
+                                                            intQuantsIn.pvtRegionIndex(),
+                                                            SoMax);
+
+        const auto& salt_in = Opm::getValue(intQuantsIn.fluidState().saltSaturation());
+
+        const auto& bLiquidIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? 
+            FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in, rs_zero_in, salt_in):
+            FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in, rs_zero_in);
+
+        const auto& bLiquidSatIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? 
+            FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in, salt_in) :
+            FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_in, p_in);
+
+        const auto& densityLiquidIn = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
+            FluidSystem::waterPvt().waterReferenceDensity(intQuantsIn.pvtRegionIndex()) :
+            FluidSystem::oilPvt().oilReferenceDensity(intQuantsIn.pvtRegionIndex());
+
+        const auto rho_in = bLiquidIn * densityLiquidIn;
+        const auto rho_sat_in = bLiquidSatIn
+            * (densityLiquidIn + rssat_in * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, intQuantsIn.pvtRegionIndex()));
+        
+        //exteriour
+        const Scalar rs_zero_ex = 0.0;
+        const auto& t_ex = Opm::getValue(intQuantsEx.fluidState().temperature(liquidPhaseIdx));
+        const auto& p_ex = Opm::getValue(intQuantsEx.fluidState().pressure(liquidPhaseIdx));
+        const auto& rssat_ex = FluidSystem::saturatedDissolutionFactor(intQuantsEx.fluidState(),
+                                                            liquidPhaseIdx,
+                                                            intQuantsEx.pvtRegionIndex(),
+                                                            SoMax);
+        const auto& salt_ex = Opm::getValue(intQuantsEx.fluidState().saltSaturation());
+            
+        const auto& bLiquidEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? 
+            FluidSystem::waterPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex, rs_zero_ex, salt_ex) :
+            FluidSystem::oilPvt().inverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex, rs_zero_ex);
+
+        const auto& bLiquidSatEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? 
+            FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex, salt_ex) :
+            FluidSystem::oilPvt().saturatedInverseFormationVolumeFactor(intQuantsIn.pvtRegionIndex(), t_ex, p_ex);
+
+        const auto& densityLiquidEx = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ? 
+            FluidSystem::waterPvt().waterReferenceDensity(intQuantsEx.pvtRegionIndex()) :
+            FluidSystem::oilPvt().oilReferenceDensity(intQuantsEx.pvtRegionIndex());
+
+
+        const auto rho_ex = bLiquidEx * densityLiquidEx;
+        const auto rho_sat_ex = bLiquidSatEx
+            * (densityLiquidEx + rssat_ex * FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, intQuantsEx.pvtRegionIndex()));
+        
+        //rho difference approximation
+        const auto delta_rho = (rho_sat_ex + rho_sat_in - rho_in -rho_ex)/2;
+        const auto pressure_difference_convective_mixing =  delta_rho * distZg;
+
+        //if change in pressure
+        if (Opm::abs(pressure_difference_convective_mixing) > 1e-12){
+
+            // find new upstream direction
+            short interiorDofIdx = 0;
+            short exteriorDofIdx = 1;
+            short upIdx = 0;
+            short downIdx = 1;
+
+            if (pressure_difference_convective_mixing > 0) {
+                upIdx = exteriorDofIdx;
+                downIdx = interiorDofIdx;
+            }
+
+            const auto& up = (upIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
+            unsigned globalUpIndex = (upIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
+
+            const auto& down = (downIdx == interiorDofIdx) ? intQuantsIn : intQuantsEx;
+            unsigned globalDownIndex = (downIdx == interiorDofIdx) ? globalIndexIn : globalIndexEx;
+
+            const auto& Rsup =  (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
+                                up.fluidState().Rsw() :
+                                up.fluidState().Rs();
+            const auto& Rsdown =  (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
+                                down.fluidState().Rsw() :
+                                down.fluidState().Rs();
+            
+            const auto& RsSat = FluidSystem::saturatedDissolutionFactor(up.fluidState(),
+                                                            liquidPhaseIdx,
+                                                            up.pvtRegionIndex(),
+                                                            SoMax);
+            const Evaluation& transMult = up.rockCompTransMultiplier();
+            
+            Evaluation sg = up.fluidState().saturation(FluidSystem::gasPhaseIdx);
+            Evaluation X = (Rsup - RsSat * sg) / (RsSat * ( 1.0 - sg));
+            if ( (X > info.Psi_[up.pvtRegionIndex()] || Rsdown > 0) ) {
+                const auto& invB = up.fluidState().invB(liquidPhaseIdx);
+                const auto& visc = up.fluidState().viscosity(liquidPhaseIdx);
+                // what will be the flux when muliplied with trans_mob
+                const auto convectiveFlux = -trans*transMult*info.Xhi_[up.pvtRegionIndex()]*invB*pressure_difference_convective_mixing*Rsup/(visc*faceArea);
+                unsigned activeGasCompIdx = Indices::canonicalToActiveComponentIndex(FluidSystem::gasCompIdx);                   
+                if (globalUpIndex == globalIndexIn)
+                    flux[conti0EqIdx + activeGasCompIdx] += convectiveFlux;
+                else
+                    flux[conti0EqIdx + activeGasCompIdx] += Opm::getValue(convectiveFlux);
+            }
+        } 
+    };
+
+};
+
+}
+#endif
+
+
+

opm/models/blackoil/blackoilintensivequantities.hh

@@ -190,6 +190,9 @@ public:
         Scalar RsMax = FluidSystem::enableDissolvedGas()
             ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
             : 0.0;
+        Scalar RswMax = FluidSystem::enableDissolvedGasInWater()
+            ? problem.maxGasDissolutionFactor(timeIdx, globalSpaceIdx)
+            : 0.0;
 
         asImp_().updateTemperature_(elemCtx, dofIdx, timeIdx);
 
@@ -300,8 +303,8 @@ public:
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
             SoMax = max(fluidState_.saturation(oilPhaseIdx),
                         problem.maxOilSaturation(globalSpaceIdx));
-        }
+        
-
+        
         // take the meaning of the switching primary variable into account for the gas
         // and oil phase compositions
         if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
@@ -310,16 +313,16 @@ public:
         } else {
             if (FluidSystem::enableDissolvedGas()) { // Add So > 0? i.e. if only water set rs = 0)
                 const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
-                    FluidSystem::saturatedDissolutionFactor(fluidState_,
+                FluidSystem::saturatedDissolutionFactor(fluidState_,
-                                                            oilPhaseIdx,
+                                                    oilPhaseIdx,
-                                                            pvtRegionIdx,
+                                                    pvtRegionIdx,
-                                                            SoMax);
+                                                    SoMax);
                 fluidState_.setRs(min(RsMax, RsSat));
             }
             else if constexpr (compositionSwitchEnabled)
                 fluidState_.setRs(0.0);
         }
-
+        }
         if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
             const auto& Rv = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
             fluidState_.setRv(Rv);
@@ -356,7 +359,7 @@ public:
                 const Evaluation& RswSat = FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             waterPhaseIdx,
                                                             pvtRegionIdx);
-                fluidState_.setRsw(RswSat);
+                fluidState_.setRsw(min(RswMax, RswSat));
             }
         }
 
@@ -402,10 +405,13 @@ public:
             rho = fluidState_.invB(waterPhaseIdx);
             rho *= FluidSystem::referenceDensity(waterPhaseIdx, pvtRegionIdx);
             if (FluidSystem::enableDissolvedGasInWater()) {
-                rho +=
+                const auto& oilVaporizationControl = problem.simulator().vanguard().schedule()[problem.episodeIndex()].oilvap();
-                    fluidState_.invB(waterPhaseIdx) *
+				if(!oilVaporizationControl.drsdtConvective()) {
-                    fluidState_.Rsw() *
+                    rho +=
-                    FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+                        fluidState_.invB(waterPhaseIdx) *
+                        fluidState_.Rsw() *
+                        FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+                }
             }
             fluidState_.setDensity(waterPhaseIdx, rho);
         }
@@ -431,11 +437,14 @@ public:
         if (FluidSystem::phaseIsActive(oilPhaseIdx)) {
             rho = fluidState_.invB(oilPhaseIdx);
             rho *= FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx);
-            if (FluidSystem::enableDissolvedGas()) {
+            if (FluidSystem::enableDissolvedGas()) {                
-                rho +=
+                const auto& oilVaporizationControl = problem.simulator().vanguard().schedule()[problem.episodeIndex()].oilvap();
-                    fluidState_.invB(oilPhaseIdx) *
+				if(!oilVaporizationControl.drsdtConvective()) {
-                    fluidState_.Rs() *
+                    rho +=
-                    FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+                        fluidState_.invB(oilPhaseIdx) *
+                        fluidState_.Rs() *
+                        FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx);
+                }
             }
             fluidState_.setDensity(oilPhaseIdx, rho);
         }

opm/models/blackoil/blackoillocalresidual.hh

@@ -37,6 +37,7 @@
 #include "blackoilbrinemodules.hh"
 #include "blackoildiffusionmodule.hh"
 #include "blackoilmicpmodules.hh"
+#include "blackoilconvectivemixingmodule.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
 
 namespace Opm {
@@ -90,6 +91,8 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
     using BrineModule = BlackOilBrineModule<TypeTag>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
     using MICPModule = BlackOilMICPModule<TypeTag>;
+    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag>;
+
 
 public:
     /*!
@@ -210,7 +213,7 @@ public:
             else
                 evalPhaseFluxes_<Scalar>(flux, phaseIdx, pvtRegionIdx, extQuants, up.fluidState());
         }
-
+		
         // deal with solvents (if present)
         SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
@@ -233,6 +236,9 @@ public:
         MICPModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
         DiffusionModule::addDiffusiveFlux(flux, elemCtx, scvfIdx, timeIdx);
+
+        // deal with convective mixing (if present)
+        ConvectiveMixingModule::addConvectiveMixingFlux(flux,elemCtx, scvfIdx, timeIdx);
     }
 
     /*!

opm/models/blackoil/blackoillocalresidualtpfa.hh

@@ -36,6 +36,7 @@
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
 #include "blackoildiffusionmodule.hh"
+#include "blackoilconvectivemixingmodule.hh"
 #include "blackoildispersionmodule.hh"
 #include "blackoilmicpmodules.hh"
 #include <opm/material/fluidstates/BlackOilFluidState.hpp>
@@ -103,6 +104,9 @@ class BlackOilLocalResidualTPFA : public GetPropType<TypeTag, Properties::DiscLo
     using FoamModule = BlackOilFoamModule<TypeTag>;
     using BrineModule = BlackOilBrineModule<TypeTag>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
+    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag>;
+    using ConvectiveMixingModuleParam = typename ConvectiveMixingModule::ConvectiveMixingModuleParam;
+
     using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
     using MICPModule = BlackOilMICPModule<TypeTag>;
 
@@ -124,6 +128,11 @@ public:
         double diffusivity;
         double dispersivity;
     };
+
+    struct ModuleParams {
+        ConvectiveMixingModuleParam convectiveMixingModuleParam;
+    };
+
     /*!
      * \copydoc FvBaseLocalResidual::computeStorage
      */
@@ -227,7 +236,8 @@ public:
                             const unsigned globalIndexEx,
                             const IntensiveQuantities& intQuantsIn,
                             const IntensiveQuantities& intQuantsEx,
-                            const ResidualNBInfo& nbInfo)
+                            const ResidualNBInfo& nbInfo,
+                            const ModuleParams& moduleParams)
     {
         OPM_TIMEBLOCK_LOCAL(computeFlux);
         flux = 0.0;
@@ -239,7 +249,8 @@ public:
                          intQuantsEx,
                          globalIndexIn,
                          globalIndexEx,
-                         nbInfo);
+                         nbInfo,
+                         moduleParams);
     }
 
     // This function demonstrates compatibility with the ElementContext-based interface.
@@ -306,7 +317,8 @@ public:
                          intQuantsEx,
                          globalIndexIn,
                          globalIndexEx,
-                         res_nbinfo);
+                         res_nbinfo,
+                         problem.moduleParams());
     }
 
     static void calculateFluxes_(RateVector& flux,
@@ -315,7 +327,8 @@ public:
                                  const IntensiveQuantities& intQuantsEx,
                                  const unsigned& globalIndexIn,
                                  const unsigned& globalIndexEx,
-                                 const ResidualNBInfo& nbInfo)
+                                 const ResidualNBInfo& nbInfo,
+                                 const ModuleParams& moduleParams)
     {
         OPM_TIMEBLOCK_LOCAL(calculateFluxes);
         const Scalar Vin = nbInfo.Vin;
@@ -395,7 +408,7 @@ public:
                 }
             }
 
-        }
+        }      
 
         // deal with solvents (if present)
         static_assert(!enableSolvent, "Relevant computeFlux() method must be implemented for this module before enabling.");
@@ -443,6 +456,18 @@ public:
         static_assert(!enableBrine, "Relevant computeFlux() method must be implemented for this module before enabling.");
         // BrineModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
+        // deal with convective mixing
+        ConvectiveMixingModule::addConvectiveMixingFlux(flux,
+                                                        intQuantsIn,
+                                                        intQuantsEx,
+                                                        globalIndexIn,
+                                                        globalIndexEx,
+                                                        nbInfo.dZg,
+                                                        nbInfo.trans,
+                                                        nbInfo.faceArea,
+                                                        moduleParams.convectiveMixingModuleParam);
+
+
         // deal with diffusion (if present). opm-models expects per area flux (added in the tmpdiffusivity).
         if constexpr(enableDiffusion){
             typename DiffusionModule::ExtensiveQuantities::EvaluationArray effectiveDiffusionCoefficient;
@@ -617,7 +642,7 @@ public:
                                        RateVector& bdyFlux,
                                        const BoundaryConditionData& bdyInfo,
                                        const IntensiveQuantities& insideIntQuants,
-                                       [[maybe_unused]] unsigned globalSpaceIdx)
+                                       unsigned globalSpaceIdx)
     {
         OPM_TIMEBLOCK_LOCAL(computeBoundaryThermal);
         // only heat is allowed to flow through this boundary

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -705,7 +705,8 @@ public:
                                                                                    pw,
                                                                                    saltConcentration);
                     setPrimaryVarsMeaningWater(WaterMeaning::Rsw);
-                    (*this)[Indices::waterSwitchIdx] = rswSat; //primary variable becomes Rsw
+                    Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
+                    (*this)[Indices::waterSwitchIdx] = min(rswSat, rswMax); //primary variable becomes Rsw
                     setPrimaryVarsMeaningPressure(PressureMeaning::Pw);
                     this->setScaledPressure_(pw);
                     changed = true;
@@ -749,7 +750,8 @@ public:
                                                                                    saltConcentration);
 
                 Scalar rsw = (*this)[Indices::waterSwitchIdx];
-                if (rsw > rswSat) {
+                Scalar rswMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
+                if (rsw > min(rswSat, rswMax)) {
                     // the gas phase appears, i.e., switch the primary variables to WaterMeaning::Sw
                     setPrimaryVarsMeaningWater(WaterMeaning::Sw);
                     (*this)[Indices::waterSwitchIdx] = 1.0; // hydrocarbon water saturation

opm/models/discretization/common/tpfalinearizer.hh

@@ -650,7 +650,7 @@ public:
             return;
         }
         const unsigned int numCells = model_().numTotalDof();
-
+        
 #ifdef _OPENMP
 #pragma omp parallel for
 #endif
@@ -671,7 +671,7 @@ public:
                 adres = 0.0;
                 darcyFlux = 0.0;
                 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
-                LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo);
+                LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo, problem_().moduleParams());
                 adres *= nbInfo.res_nbinfo.faceArea;
                 if (enableFlows) {
                     for (unsigned eqIdx = 0; eqIdx < numEq; ++ eqIdx) {
@@ -757,7 +757,7 @@ private:
                 adres = 0.0;
                 darcyFlux = 0.0;
                 const IntensiveQuantities& intQuantsEx = model_().intensiveQuantities(globJ, /*timeIdx*/ 0);
-                LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo);
+                LocalResidual::computeFlux(adres,darcyFlux, globI, globJ, intQuantsIn, intQuantsEx, nbInfo.res_nbinfo,  problem_().moduleParams());
                 adres *= nbInfo.res_nbinfo.faceArea;
                 if (enableDispersion) {
                     for (unsigned phaseIdx = 0; phaseIdx < numEq; ++ phaseIdx) {