Alternative solvent extension for the black oil model.

opm/models/blackoil/blackoilindices.hh

@@ -35,7 +35,7 @@ namespace Opm {
  *
  * \brief The primary variable and equation indices for the black-oil model.
  */
-template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset>
+template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset>
 struct BlackOilIndices
 {
     //! Number of phases active at all times
@@ -49,6 +49,9 @@ struct BlackOilIndices
     //! Are solvents involved?
     static const bool enableSolvent = numSolventsV > 0;
 
+    //! Is extbo invoked?
+    static const bool enableExtbo = numExtbosV > 0;
+
     //! Are polymers involved?
     static const bool enablePolymer = numPolymersV > 0;
 
@@ -58,6 +61,9 @@ struct BlackOilIndices
     //! Number of solvent components to be considered
     static const int numSolvents = enableSolvent ? numSolventsV : 0;
 
+    //! Number of components to be considered for extbo
+    static const int numExtbos = enableExtbo ? numExtbosV : 0;
+
     //! Number of polymer components to be considered
     static const int numPolymers = enablePolymer ? numPolymersV : 0;
 
@@ -71,7 +77,7 @@ struct BlackOilIndices
     static const int numBrine = enableBrine? 1 : 0;
 
     //! The number of equations
-    static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam + numBrine;
+    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
 
     //! \brief returns the index of "active" component
     static constexpr unsigned canonicalToActiveComponentIndex(unsigned compIdx)
@@ -104,6 +110,10 @@ struct BlackOilIndices
     static const int solventSaturationIdx =
         enableSolvent ? PVOffset + numPhases : -1000;
 
+    //! Index of the primary variable for the first extbo component
+    static const int zFractionIdx =
+        enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
+
     //! Index of the primary variable for the first polymer
     static const int polymerConcentrationIdx =
         enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@@ -114,15 +124,15 @@ struct BlackOilIndices
 
     //! Index of the primary variable for the foam
     static const int foamConcentrationIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000;
 
     //! Index of the primary variable for the brine
     static const int saltConcentrationIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numExtbos + numPolymers + numFoam : -1000;
 
     //! Index of the primary variable for temperature
     static const int temperatureIdx  =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine : - 1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
 
 
     ////////
@@ -137,9 +147,13 @@ struct BlackOilIndices
     static const int contiSolventEqIdx =
         enableSolvent ? PVOffset + numPhases : -1000;
 
+    //! Index of the continuity equation for the first extbo component
+    static const int contiZfracEqIdx =
+        enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
+
     //! Index of the continuity equation for the first polymer component
     static const int contiPolymerEqIdx =
-        enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
+        enablePolymer ? PVOffset + numPhases + numSolvents + numExtbos : -1000;
 
     //! Index of the continuity equation for the second polymer component (molecular weight)
     static const int contiPolymerMWEqIdx =
@@ -147,16 +161,16 @@ struct BlackOilIndices
 
     //! Index of the continuity equation for the foam component
     static const int contiFoamEqIdx =
-        enableFoam ? PVOffset + numPhases + numSolvents + numPolymers : -1000;
+        enableFoam ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers : -1000;
 
     //! Index of the continuity equation for the salt water component
     static const int contiBrineEqIdx =
-        enableBrine ? PVOffset + numPhases + numSolvents + numPolymers + numFoam : -1000;
+        enableBrine ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam : -1000;
 
 
     //! Index of the continuity equation for energy
     static const int contiEnergyEqIdx =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine: -1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine: -1000;
 };
 
 } // namespace Opm

opm/models/blackoil/blackoilintensivequantities.hh

@@ -30,6 +30,7 @@
 
 #include "blackoilproperties.hh"
 #include "blackoilsolventmodules.hh"
+#include "blackoilextbomodules.hh"
 #include "blackoilpolymermodules.hh"
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
@@ -55,6 +56,7 @@ class BlackOilIntensiveQuantities
     : public GetPropType<TypeTag, Properties::DiscIntensiveQuantities>
     , public GetPropType<TypeTag, Properties::FluxModule>::FluxIntensiveQuantities
     , public BlackOilSolventIntensiveQuantities<TypeTag>
+    , public BlackOilExtboIntensiveQuantities<TypeTag>
     , public BlackOilPolymerIntensiveQuantities<TypeTag>
     , public BlackOilFoamIntensiveQuantities<TypeTag>
     , public BlackOilBrineIntensiveQuantities<TypeTag>
@@ -75,6 +77,7 @@ class BlackOilIntensiveQuantities
 
     enum { numEq = getPropValue<TypeTag, Properties::NumEq>() };
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+    enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
     enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
@@ -210,6 +213,9 @@ public:
         // update the Saturation functions for the blackoil solvent module.
         asImp_().solventPostSatFuncUpdate_(elemCtx, dofIdx, timeIdx);
 
+        // update extBO parameters
+        asImp_().zFractionUpdate_(elemCtx, dofIdx, timeIdx);
+
         Evaluation SoMax = 0.0;
         if (FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
             SoMax = Opm::max(fluidState_.saturation(oilPhaseIdx),
@@ -223,7 +229,7 @@ public:
             // we use the compositions of the gas-saturated oil and oil-saturated gas.
             if (FluidSystem::enableDissolvedGas()) {
                 Scalar RsMax = elemCtx.problem().maxGasDissolutionFactor(timeIdx, globalSpaceIdx);
-                const Evaluation& RsSat =
+                const Evaluation& RsSat = enableExtbo ? asImp_().rs() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             oilPhaseIdx,
                                                             pvtRegionIdx,
@@ -235,7 +241,7 @@ public:
 
             if (FluidSystem::enableVaporizedOil()) {
                 Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
-                const Evaluation& RvSat =
+                const Evaluation& RvSat = enableExtbo ? asImp_().rv() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             gasPhaseIdx,
                                                             pvtRegionIdx,
@@ -257,7 +263,7 @@ public:
                 // the gas phase is not present, but we need to compute its "composition"
                 // for the gravity correction anyway
                 Scalar RvMax = elemCtx.problem().maxOilVaporizationFactor(timeIdx, globalSpaceIdx);
-                const auto& RvSat =
+                const auto& RvSat = enableExtbo ? asImp_().rv() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             gasPhaseIdx,
                                                             pvtRegionIdx,
@@ -276,7 +282,7 @@ public:
                 // 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 =
+                const auto& RsSat = enableExtbo ? asImp_().rs() :
                     FluidSystem::saturatedDissolutionFactor(fluidState_,
                                                             oilPhaseIdx,
                                                             pvtRegionIdx,
@@ -306,7 +312,12 @@ public:
             fluidState_.setInvB(phaseIdx, b);
 
             const auto& mu = FluidSystem::viscosity(fluidState_, paramCache, phaseIdx);
-            mobility_[phaseIdx] /= mu;
+            if (enableExtbo && phaseIdx == oilPhaseIdx)
+              mobility_[phaseIdx] /= asImp_().visco();
+            else if (enableExtbo && phaseIdx == gasPhaseIdx)
+              mobility_[phaseIdx] /= asImp_().viscg();
+            else
+              mobility_[phaseIdx] /= mu;
         }
         Opm::Valgrind::CheckDefined(mobility_);
 
@@ -366,6 +377,7 @@ public:
         porosity_ *= problem.template rockCompPoroMultiplier<Evaluation>(*this, globalSpaceIdx);
 
         asImp_().solventPvtUpdate_(elemCtx, dofIdx, timeIdx);
+        asImp_().zPvtUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().polymerPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
         asImp_().updateEnergyQuantities_(elemCtx, dofIdx, timeIdx, paramCache);
         asImp_().foamPropertiesUpdate_(elemCtx, dofIdx, timeIdx);
@@ -446,6 +458,7 @@ public:
 
 private:
     friend BlackOilSolventIntensiveQuantities<TypeTag>;
+    friend BlackOilExtboIntensiveQuantities<TypeTag>;
     friend BlackOilPolymerIntensiveQuantities<TypeTag>;
     friend BlackOilEnergyIntensiveQuantities<TypeTag>;
     friend BlackOilFoamIntensiveQuantities<TypeTag>;

opm/models/blackoil/blackoillocalresidual.hh

@@ -30,6 +30,7 @@
 
 #include "blackoilproperties.hh"
 #include "blackoilsolventmodules.hh"
+#include "blackoilextbomodules.hh"
 #include "blackoilpolymermodules.hh"
 #include "blackoilenergymodules.hh"
 #include "blackoilfoammodules.hh"
@@ -80,6 +81,7 @@ class BlackOilLocalResidual : public GetPropType<TypeTag, Properties::DiscLocalR
 
     using Toolbox = Opm::MathToolbox<Evaluation>;
     using SolventModule = BlackOilSolventModule<TypeTag>;
+    using ExtboModule = BlackOilExtboModule<TypeTag>;
     using PolymerModule = BlackOilPolymerModule<TypeTag>;
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
     using FoamModule = BlackOilFoamModule<TypeTag>;
@@ -143,6 +145,9 @@ public:
         // deal with solvents (if present)
         SolventModule::addStorage(storage, intQuants);
 
+        // deal with zFracton (if present)
+        ExtboModule::addStorage(storage, intQuants);
+
         // deal with polymer (if present)
         PolymerModule::addStorage(storage, intQuants);
 
@@ -186,6 +191,9 @@ public:
         // deal with solvents (if present)
         SolventModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 
+        // deal with zFracton (if present)
+        ExtboModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
+
         // deal with polymer (if present)
         PolymerModule::computeFlux(flux, elemCtx, scvfIdx, timeIdx);
 

opm/models/blackoil/blackoilmodel.hh

@@ -45,6 +45,7 @@
 #include "blackoilpolymermodules.hh"
 #include "blackoilfoammodules.hh"
 #include "blackoilbrinemodules.hh"
+#include "blackoilextbomodules.hh"
 #include "blackoildarcyfluxmodule.hh"
 
 #include <opm/models/common/multiphasebasemodel.hh>
@@ -122,6 +123,7 @@ struct FluxModule<TypeTag, TTag::BlackOilModel> { using type = Opm::BlackOilDarc
 template<class TypeTag>
 struct Indices<TypeTag, TTag::BlackOilModel>
 { using type = Opm::BlackOilIndices<getPropValue<TypeTag, Properties::EnableSolvent>(),
+                                    getPropValue<TypeTag, Properties::EnableExtbo>(),
                                     getPropValue<TypeTag, Properties::EnablePolymer>(),
                                     getPropValue<TypeTag, Properties::EnableEnergy>(),
                                     getPropValue<TypeTag, Properties::EnableFoam>(),
@@ -144,6 +146,8 @@ public:
 template<class TypeTag>
 struct EnableSolvent<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
 template<class TypeTag>
+struct EnableExtbo<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
+template<class TypeTag>
 struct EnablePolymer<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
 template<class TypeTag>
 struct EnablePolymerMW<TypeTag, TTag::BlackOilModel> { static constexpr bool value = false; };
@@ -271,6 +275,7 @@ class BlackOilModel
     static const bool waterEnabled = Indices::waterEnabled;
 
     using SolventModule = BlackOilSolventModule<TypeTag>;
+    using ExtboModule = BlackOilExtboModule<TypeTag>;
     using PolymerModule = BlackOilPolymerModule<TypeTag>;
     using EnergyModule = BlackOilEnergyModule<TypeTag>;
 public:
@@ -286,6 +291,7 @@ public:
         ParentType::registerParameters();
 
         SolventModule::registerParameters();
+        ExtboModule::registerParameters();
         PolymerModule::registerParameters();
         EnergyModule::registerParameters();
 
@@ -315,6 +321,8 @@ public:
             oss << "composition_switching";
         else if (SolventModule::primaryVarApplies(pvIdx))
             return SolventModule::primaryVarName(pvIdx);
+        else if (ExtboModule::primaryVarApplies(pvIdx))
+            return ExtboModule::primaryVarName(pvIdx);
         else if (PolymerModule::primaryVarApplies(pvIdx))
             return PolymerModule::primaryVarName(pvIdx);
         else if (EnergyModule::primaryVarApplies(pvIdx))
@@ -336,6 +344,8 @@ public:
             oss << "conti_" << FluidSystem::phaseName(eqIdx - Indices::conti0EqIdx);
         else if (SolventModule::eqApplies(eqIdx))
             return SolventModule::eqName(eqIdx);
+        else if (ExtboModule::eqApplies(eqIdx))
+            return ExtboModule::eqName(eqIdx);
         else if (PolymerModule::eqApplies(eqIdx))
             return PolymerModule::eqName(eqIdx);
         else if (EnergyModule::eqApplies(eqIdx))
@@ -369,6 +379,10 @@ public:
         else if (SolventModule::primaryVarApplies(pvIdx))
             return SolventModule::primaryVarWeight(pvIdx);
 
+        // deal with primary variables stemming from the extBO module
+        else if (ExtboModule::primaryVarApplies(pvIdx))
+            return ExtboModule::primaryVarWeight(pvIdx);
+
         // deal with primary variables stemming from the polymer module
         else if (PolymerModule::primaryVarApplies(pvIdx))
             return PolymerModule::primaryVarWeight(pvIdx);
@@ -424,6 +438,9 @@ public:
         if (SolventModule::eqApplies(eqIdx))
             return SolventModule::eqWeight(eqIdx);
 
+        if (ExtboModule::eqApplies(eqIdx))
+            return ExtboModule::eqWeight(eqIdx);
+
         else if (PolymerModule::eqApplies(eqIdx))
             return PolymerModule::eqWeight(eqIdx);
 
@@ -463,6 +480,7 @@ public:
         outstream << priVars.pvtRegionIndex() << " ";
 
         SolventModule::serializeEntity(*this, outstream, dof);
+        ExtboModule::serializeEntity(*this, outstream, dof);
         PolymerModule::serializeEntity(*this, outstream, dof);
         EnergyModule::serializeEntity(*this, outstream, dof);
     }
@@ -500,6 +518,7 @@ public:
             throw std::runtime_error("Could not deserialize degree of freedom "+std::to_string(dofIdx));
 
         SolventModule::deserializeEntity(*this, instream, dof);
+        ExtboModule::deserializeEntity(*this, instream, dof);
         PolymerModule::deserializeEntity(*this, instream, dof);
         EnergyModule::deserializeEntity(*this, instream, dof);
 

opm/models/blackoil/blackoilnewtonmethod.hh

@@ -212,6 +212,7 @@ protected:
                                  const EqVector& currentResidual)
     {
         static constexpr bool enableSolvent = Indices::solventSaturationIdx >= 0;
+        static constexpr bool enableExtbo = Indices::zFractionIdx >= 0;
         static constexpr bool enablePolymer = Indices::polymerConcentrationIdx >= 0;
         static constexpr bool enablePolymerWeight = Indices::polymerMoleWeightIdx >= 0;
         static constexpr bool enableEnergy = Indices::temperatureIdx >= 0;
@@ -287,6 +288,13 @@ protected:
             else if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
                 // solvent saturation updates are also subject to the Appleyard chop
                 delta *= satAlpha;
+            else if (enableExtbo && pvIdx == Indices::zFractionIdx) {
+                // z fraction updates are also subject to the Appleyard chop
+                if (delta > currentValue[Indices::zFractionIdx])
+                        delta = currentValue[Indices::zFractionIdx];
+                if (delta < currentValue[Indices::zFractionIdx]-1.0)
+                        delta = currentValue[Indices::zFractionIdx]-1.0;
+            }
             else if (enablePolymerWeight && pvIdx == Indices::polymerMoleWeightIdx) {
                 const double sign = delta >= 0. ? 1. : -1.;
                 // maximum change of polymer molecular weight, the unit is MDa.
@@ -303,6 +311,10 @@ protected:
             if (enableSolvent && pvIdx == Indices::solventSaturationIdx)
                 nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
 
+            // keep the z fraction between 0 and 1
+            if (enableExtbo && pvIdx == Indices::zFractionIdx)
+                nextValue[pvIdx] = std::min(std::max(nextValue[pvIdx], 0.0), 1.0);
+
             // keep the polymer concentration above 0
             if (enablePolymer && pvIdx == Indices::polymerConcentrationIdx)
                 nextValue[pvIdx] = std::max(nextValue[pvIdx], 0.0);

opm/models/blackoil/blackoilprimaryvariables.hh

@@ -30,6 +30,7 @@
 
 #include "blackoilproperties.hh"
 #include "blackoilsolventmodules.hh"
+#include "blackoilextbomodules.hh"
 #include "blackoilpolymermodules.hh"
 #include "blackoilenergymodules.hh"
 #include "blackoilfoammodules.hh"
@@ -49,6 +50,9 @@ namespace Opm {
 template <class TypeTag, bool enableSolvent>
 class BlackOilSolventModule;
 
+template <class TypeTag, bool enableExtbo>
+class BlackOilExtboModule;
+
 template <class TypeTag, bool enablePolymer>
 class BlackOilPolymerModule;
 
@@ -94,6 +98,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     // component indices from the fluid system
     enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
     enum { enableSolvent = getPropValue<TypeTag, Properties::EnableSolvent>() };
+    enum { enableExtbo = getPropValue<TypeTag, Properties::EnableExtbo>() };
     enum { enablePolymer = getPropValue<TypeTag, Properties::EnablePolymer>() };
     enum { enableFoam = getPropValue<TypeTag, Properties::EnableFoam>() };
     enum { enableBrine = getPropValue<TypeTag, Properties::EnableBrine>() };
@@ -105,6 +110,7 @@ class BlackOilPrimaryVariables : public FvBasePrimaryVariables<TypeTag>
     using Toolbox = typename Opm::MathToolbox<Evaluation>;
     using ComponentVector = Dune::FieldVector<Scalar, numComponents>;
     using SolventModule = BlackOilSolventModule<TypeTag, enableSolvent>;
+    using ExtboModule = BlackOilExtboModule<TypeTag, enableExtbo>;
     using PolymerModule = BlackOilPolymerModule<TypeTag, enablePolymer>;
     using EnergyModule = BlackOilEnergyModule<TypeTag, enableEnergy>;
     using FoamModule = BlackOilFoamModule<TypeTag, enableFoam>;
@@ -402,7 +408,10 @@ public:
                 Scalar T = asImp_().temperature_();
                 Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
                 Scalar RsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
-                Scalar RsSat = FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
+                Scalar RsSat = enableExtbo ? ExtboModule::rs(pvtRegionIndex(),
+                                                             po,
+                                                             zFraction_())
+                             : FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
                                                                                    T,
                                                                                    po,
                                                                                    So2,
@@ -412,6 +421,7 @@ public:
                 if (compositionSwitchEnabled)
                     (*this)[Indices::compositionSwitchIdx] =
                         std::min(RsMax, RsSat);
+
                 return true;
             }
 
@@ -433,12 +443,14 @@ public:
                 Scalar T = asImp_().temperature_();
                 Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
                 Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
-                Scalar RvSat =
+                Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
-                    FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                             pg,
-                                                                         T,
+                                                             zFraction_())
-                                                                         pg,
+                             : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
-                                                                         Scalar(0),
+                                                                                    T,
-                                                                         SoMax);
+                                                                                    pg,
+                                                                                    Scalar(0),
+                                                                                    SoMax);
                 setPrimaryVarsMeaning(Sw_pg_Rv);
                 (*this)[Indices::pressureSwitchIdx] = pg;
                 if (compositionSwitchEnabled)
@@ -473,12 +485,14 @@ public:
             Scalar So = 1.0 - Sw - solventSaturation_();
             Scalar SoMax = std::max(So, problem.maxOilSaturation(globalDofIdx));
             Scalar RsMax = problem.maxGasDissolutionFactor(/*timeIdx=*/0, globalDofIdx);
-            Scalar RsSat =
+            Scalar RsSat = enableExtbo ? ExtboModule::rs(pvtRegionIndex(),
-                FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
+                                                         po,
-                                                                    T,
+                                                         zFraction_())
-                                                                    po,
+                         : FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_,
-                                                                    So,
+                                                                               T,
-                                                                    SoMax);
+                                                                               po,
+                                                                               So,
+                                                                               SoMax);
 
             Scalar Rs = (*this)[Indices::compositionSwitchIdx];
             if (Rs > std::min(RsMax, RsSat*(1.0 + eps))) {
@@ -530,12 +544,14 @@ public:
             Scalar T = asImp_().temperature_();
             Scalar SoMax = problem.maxOilSaturation(globalDofIdx);
             Scalar RvMax = problem.maxOilVaporizationFactor(/*timeIdx=*/0, globalDofIdx);
-            Scalar RvSat =
+            Scalar RvSat = enableExtbo ? ExtboModule::rv(pvtRegionIndex(),
-                FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
+                                                         pg,
-                                                                     T,
+                                                         zFraction_())
-                                                                     pg,
+                         : FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_,
-                                                                     /*So=*/Scalar(0.0),
+                                                                                T,
-                                                                     SoMax);
+                                                                                pg,
+                                                                                /*So=*/Scalar(0.0),
+                                                                                SoMax);
 
             Scalar Rv = (*this)[Indices::compositionSwitchIdx];
             if (Rv > std::min(RvMax, RvSat*(1.0 + eps))) {
@@ -686,6 +702,14 @@ private:
         return (*this)[Indices::solventSaturationIdx];
     }
 
+    Scalar zFraction_() const
+    {
+        if (!enableExtbo)
+            return 0.0;
+
+        return (*this)[Indices::zFractionIdx];
+    }
+
     Scalar polymerConcentration_() const
     {
         if (!enablePolymer)

opm/models/blackoil/blackoilproperties.hh

@@ -40,6 +40,9 @@ struct EnableEclipseOutput { using type = UndefinedProperty; };
 //! Enable the ECL-blackoil extension for solvents. ("Second gas")
 template<class TypeTag, class MyTypeTag>
 struct EnableSolvent { using type = UndefinedProperty; };
+//! Enable the ECL-blackoil extension for extended BO. ("Second gas" - alternative approach)
+template<class TypeTag, class MyTypeTag>
+struct EnableExtbo { using type = UndefinedProperty; };
 //! Enable the ECL-blackoil extension for polymer.
 template<class TypeTag, class MyTypeTag>
 struct EnablePolymer { using type = UndefinedProperty; };

opm/models/blackoil/blackoiltwophaseindices.hh

@@ -37,7 +37,7 @@ namespace Opm {
  *
  * \brief The primary variable and equation indices for the black-oil model.
  */
-template <unsigned numSolventsV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned disabledCanonicalCompIdx>
+template <unsigned numSolventsV, unsigned numExtbosV, unsigned numPolymersV, unsigned numEnergyV, bool enableFoam, bool enableBrine, unsigned PVOffset, unsigned disabledCanonicalCompIdx>
 struct BlackOilTwoPhaseIndices
 {
     //! Is phase enabled or not
@@ -48,6 +48,9 @@ struct BlackOilTwoPhaseIndices
     //! Are solvents involved?
     static const bool enableSolvent = numSolventsV > 0;
 
+    //! Is extbo invoked?
+    static const bool enableExtbo = numExtbosV > 0;
+
     //! Are polymers involved?
     static const bool enablePolymer = numPolymersV > 0;
 
@@ -57,6 +60,9 @@ struct BlackOilTwoPhaseIndices
     //! Number of solvent components to be considered
     static const int numSolvents = enableSolvent ? numSolventsV : 0;
 
+    //! Number of components to be considered for extbo
+    static const int numExtbos = enableExtbo ? numExtbosV : 0;
+
     //! Number of polymer components to be considered
     static const int numPolymers = enablePolymer ? numPolymersV : 0;
 
@@ -73,7 +79,7 @@ struct BlackOilTwoPhaseIndices
     static const int numPhases = 2;
 
     //! The number of equations
-    static const int numEq = numPhases + numSolvents + numPolymers + numEnergy + numFoam + numBrine;
+    static const int numEq = numPhases + numSolvents + numExtbos + numPolymers + numEnergy + numFoam + numBrine;
 
     //////////////////////////////
     // Primary variable indices
@@ -97,6 +103,10 @@ struct BlackOilTwoPhaseIndices
     static const int solventSaturationIdx =
         enableSolvent ? PVOffset + numPhases : -1000;
 
+    //! Index of the primary variable for the first extbo component
+    static const int zFractionIdx =
+        enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
+
     //! Index of the primary variable for the first polymer
     static const int polymerConcentrationIdx =
         enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@@ -115,7 +125,7 @@ struct BlackOilTwoPhaseIndices
 
     //! Index of the primary variable for temperature
     static const int temperatureIdx  =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine : - 1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : - 1000;
 
     //////////////////////
     // Equation indices
@@ -166,6 +176,10 @@ struct BlackOilTwoPhaseIndices
     static const int contiSolventEqIdx =
         enableSolvent ? PVOffset + numPhases : -1000;
 
+    //! Index of the continuity equation for the first extbo component
+    static const int contiZfracEqIdx =
+        enableExtbo ? PVOffset + numPhases + numSolvents : -1000;
+
     //! Index of the continuity equation for the first polymer component
     static const int contiPolymerEqIdx =
         enablePolymer ? PVOffset + numPhases + numSolvents : -1000;
@@ -184,7 +198,7 @@ struct BlackOilTwoPhaseIndices
 
     //! Index of the continuity equation for energy
     static const int contiEnergyEqIdx =
-        enableEnergy ? PVOffset + numPhases + numSolvents + numPolymers + numFoam + numBrine : -1000;
+        enableEnergy ? PVOffset + numPhases + numSolvents + numExtbos + numPolymers + numFoam + numBrine : -1000;
 };
 
 } // namespace Opm