adressed minor comments by reviewer, added enableConvectiveMixing flag to drsdtcon, fixed indentation and bug in lastRs for gas/water

opm/simulators/flow/BlackoilModel.hpp

@@ -127,8 +127,11 @@ struct EnableMICP<TypeTag, TTag::FlowProblem>
 { static constexpr bool value = false; };
 
 template<class TypeTag>
-struct EnableDispersion<TypeTag, TTag::FlowProblem>
+struct EnableDispersion<TypeTag, TTag::FlowProblem> 
 { static constexpr bool value = false; };
+template<class TypeTag>
+struct EnableConvectiveMixing<TypeTag, TTag::FlowProblem> 
+{ static constexpr bool value = true; };
 
 template<class TypeTag>
 struct WellModel<TypeTag, TTag::FlowProblem>

opm/simulators/flow/EquilInitializer.hpp

@@ -174,7 +174,6 @@ public:
                     const auto& h = FluidSystem::enthalpy(fluidState, phaseIdx, regionIdx);
                     fluidState.setEnthalpy(phaseIdx, h);
                 }
-
             }
 
             // set the salt concentration

opm/simulators/flow/FlowProblem.hpp

@@ -145,6 +145,7 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     enum { enableEnergy = getPropValue<TypeTag, Properties::EnableEnergy>() };
     enum { enableDiffusion = getPropValue<TypeTag, Properties::EnableDiffusion>() };
     enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
+    enum { enableConvectiveMixing = getPropValue<TypeTag, Properties::EnableConvectiveMixing>() };
     enum { enableThermalFluxBoundaries = getPropValue<TypeTag, Properties::EnableThermalFluxBoundaries>() };
     enum { enableApiTracking = getPropValue<TypeTag, Properties::EnableApiTracking>() };
     enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
@@ -182,7 +183,7 @@ class FlowProblem : public GetPropType<TypeTag, Properties::BaseProblem>
     using MICPModule = BlackOilMICPModule<TypeTag>;
     using DispersionModule = BlackOilDispersionModule<TypeTag, enableDispersion>;
     using DiffusionModule = BlackOilDiffusionModule<TypeTag, enableDiffusion>;
-    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag>;
+    using ConvectiveMixingModule = BlackOilConvectiveMixingModule<TypeTag, enableConvectiveMixing>;
     using ModuleParams = typename BlackOilLocalResidualTPFA<TypeTag>::ModuleParams;
 
     using InitialFluidState = typename EquilInitializer<TypeTag>::ScalarFluidState;
@@ -538,8 +539,6 @@ public:
         const auto& schedule = simulator.vanguard().schedule();
         const auto& events = schedule[episodeIdx].events();
 
-
-
         if (episodeIdx >= 0 && events.hasEvent(ScheduleEvents::GEO_MODIFIER)) {
             // bring the contents of the keywords to the current state of the SCHEDULE
             // section.

opm/simulators/flow/FlowProblemProperties.hpp

@@ -198,7 +198,12 @@ template<class TypeTag>
 struct EnableDispersion<TypeTag, TTag::FlowBaseProblem>
 { static constexpr bool value = false; };
 
-// disable API tracking
+// Enable Convective Mixing
+template<class TypeTag>
+struct EnableConvectiveMixing<TypeTag, TTag::FlowBaseProblem>
+{ static constexpr bool value = true; };
+
+// only write the solutions for the report steps to disk
 template<class TypeTag>
 struct EnableApiTracking<TypeTag, TTag::FlowBaseProblem>
 { static constexpr bool value = false; };

opm/simulators/flow/MixingRateControls.cpp

@@ -108,7 +108,6 @@ init(std::size_t numDof, int episodeIdx, const unsigned ntpvt)
         dRsDtOnlyFreeGas_.resize(ntpvt, false);
         lastRs_.resize(numDof, 0.0);
         maxDRv_.resize(ntpvt, 1e30);
-        lastRv_.resize(numDof, 0.0);
         if (this->drsdtConvective(episodeIdx)) {
             convectiveDrs_.resize(numDof, 1.0);
         }
@@ -119,8 +118,6 @@ bool MixingRateControls<FluidSystem>::
 drsdtActive(int episodeIdx) const
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
     return (oilVaporizationControl.drsdtActive());
 }
 
@@ -129,8 +126,6 @@ bool MixingRateControls<FluidSystem>::
 drvdtActive(int episodeIdx) const
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
     return (oilVaporizationControl.drvdtActive());
 }
 
@@ -139,8 +134,6 @@ bool MixingRateControls<FluidSystem>::
 drsdtConvective(int episodeIdx) const
 {
     const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-    const bool bothOilGasActive = FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
-                                  FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx);
     return (oilVaporizationControl.drsdtConvective());
 }
 
@@ -269,8 +262,8 @@ void MixingRateControls<FluidSystem>::
 updateConvectiveDRsDt_(const unsigned compressedDofIdx,
                        const Scalar t,
                        const Scalar p,
+                       const Scalar pg,
                        const Scalar rs,
-                       const Scalar so,
                        const Scalar sg,
                        const Scalar poro,
                        const Scalar permz,
@@ -283,7 +276,7 @@ updateConvectiveDRsDt_(const unsigned compressedDofIdx,
                        const int pvtRegionIndex)
 {
     const Scalar rssat = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
-        FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p, salt) : 
+        FluidSystem::waterPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p, salt) :
         FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIndex, t, p);
     const Scalar saturatedInvB = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
         FluidSystem::waterPvt().saturatedInverseFormationVolumeFactor(pvtRegionIndex, t, p, salt) :
@@ -305,7 +298,7 @@ updateConvectiveDRsDt_(const unsigned compressedDofIdx,
         FluidSystem::waterPvt().viscosity(pvtRegionIndex, t, p, rs, salt) :
         FluidSystem::oilPvt().viscosity(pvtRegionIndex, t, p, rs);
 
-    // Note that for so = 0 this gives no limits (inf) for the dissolution rate
+    // Note that for sLiquid = 0 this gives no limits (inf) for the dissolution rate
     // Also we restrict the effect of convective mixing to positive density differences
     // i.e. we only allow for fingers moving downward
 
@@ -314,18 +307,19 @@ updateConvectiveDRsDt_(const unsigned compressedDofIdx,
 	Scalar factor = 1.0;
 	Scalar X = (rs - rssat * sg) / (rssat * ( 1.0 - sg));
     Scalar omega = 0.0;
-	if ((rs >= (rssat * sg))){
+    const Scalar pCap = Opm::abs(pg - p);
+	if ((rs >= (rssat * sg)) || (pCap < 1e-12)){
 	    if(X > Psi){
 		    factor = 0.0;
             omega = omegainn;
         }
 	} else {
-	        factor /= Xhi;
+	    factor /= Xhi;
-	        deltaDensity = (saturatedDensity - co2Density);
+	    deltaDensity = (saturatedDensity - co2Density);
-	    }
+	}
     
     convectiveDrs_[compressedDofIdx]
-        = factor * permz * rssat * max(0.0, deltaDensity) * gravity / ( std::max(sg_max - sg, 0.0) * visc * distZ * poro) + (omega/Xhi); 
+        = factor * permz * rssat * max(0.0, deltaDensity) * gravity / ( std::max(sg_max - sg, 0.0) * visc * distZ * poro) + (omega/Xhi);
 }
 
 template class MixingRateControls<BlackOilFluidSystem<double,BlackOilDefaultIndexTraits>>;

opm/simulators/flow/MixingRateControls.hpp

@@ -118,7 +118,6 @@ public:
                 const std::array<bool,3>& active)
     {
         const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
-
         if (active[0]) {
             // This implements the convective DRSDT as described in
             // Sandve et al. "Convective dissolution in field scale CO2 storage simulations using the OPM Flow
@@ -132,11 +131,7 @@ public:
             const auto& pressure = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
                 getValue(fs.pressure(FluidSystem::waterPhaseIdx)) :
                 getValue(fs.pressure(FluidSystem::oilPhaseIdx));
-
+            const auto& pressuregas = getValue(fs.pressure(FluidSystem::gasPhaseIdx));
-            const auto& sLiquid = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
-                getValue(fs.saturation(FluidSystem::waterPhaseIdx)) :
-                getValue(fs.saturation(FluidSystem::oilPhaseIdx));
-            
             const auto& rs = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
                 getValue(fs.Rsw()) :
                 getValue(fs.Rs());
@@ -146,8 +141,8 @@ public:
             this->updateConvectiveDRsDt_(compressedDofIdx,
                                          temperature,
                                          pressure,
+                                         pressuregas,
                                          rs,
-                                         sLiquid,
                                          getValue(fs.saturation(FluidSystem::gasPhaseIdx)),
                                          getValue(iq.porosity()),
                                          permZ,
@@ -164,12 +159,11 @@ public:
             const auto& fs = iq.fluidState();
 
             using FluidState = typename std::decay<decltype(fs)>::type;
-            const auto& oilVaporizationControl = schedule_[episodeIdx].oilvap();
             constexpr Scalar freeGasMinSaturation_ = 1e-7;
             if (oilVaporizationControl.getOption(pvtRegionIdx) ||
                 fs.saturation(FluidSystem::gasPhaseIdx) > freeGasMinSaturation_) {
                 lastRs_[compressedDofIdx]
-                    = (FluidSystem::phaseIsActive(FluidSystem::waterPhaseIdx)) ?
+                    = ((FluidSystem::enableDissolvedGasInWater())) ?
                     BlackOil::template getRsw_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex()) :
                     BlackOil::template getRs_<FluidSystem, FluidState, Scalar>(fs, iq.pvtRegionIndex());
             }
@@ -189,8 +183,8 @@ private:
     void updateConvectiveDRsDt_(const unsigned compressedDofIdx,
                                 const Scalar t,
                                 const Scalar p,
+                                const Scalar pg,
                                 const Scalar rs,
-                                const Scalar so,
                                 const Scalar sg,
                                 const Scalar poro,
                                 const Scalar permz,