BlackOilFluidSystem.hpp

Go to the documentation of this file.

// -*- 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.
*/
#ifndef OPM_BLACK_OIL_FLUID_SYSTEM_HPP
#define OPM_BLACK_OIL_FLUID_SYSTEM_HPP
 
#include "BlackOilDefaultIndexTraits.hpp"
#include "blackoilpvt/OilPvtMultiplexer.hpp"
#include "blackoilpvt/GasPvtMultiplexer.hpp"
#include "blackoilpvt/WaterPvtMultiplexer.hpp"
#include "blackoilpvt/BrineCo2Pvt.hpp"
 
#include <opm/material/fluidsystems/BaseFluidSystem.hpp>
#include <opm/material/Constants.hpp>
 
#include <opm/material/common/MathToolbox.hpp>
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/HasMemberGeneratorMacros.hpp>
#include <opm/material/common/Exceptions.hpp>
 
#if HAVE_ECL_INPUT
#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
#include <opm/input/eclipse/EclipseState/Tables/FlatTable.hpp>
#include <opm/input/eclipse/EclipseState/Tables/TableManager.hpp>
#endif
 
#include <memory>
#include <vector>
#include <array>
 
namespace Opm {
namespace BlackOil {
OPM_GENERATE_HAS_MEMBER(Rs, ) // Creates 'HasMember_Rs<T>'.
OPM_GENERATE_HAS_MEMBER(Rv, ) // Creates 'HasMember_Rv<T>'.
OPM_GENERATE_HAS_MEMBER(Rvw, ) // Creates 'HasMember_Rvw<T>'.
OPM_GENERATE_HAS_MEMBER(saltConcentration, )
OPM_GENERATE_HAS_MEMBER(saltSaturation, )
 
template <class FluidSystem, class FluidState, class LhsEval>
LhsEval getRs_(typename std::enable_if<!HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
               unsigned regionIdx)
{
    const auto& XoG =
        decay<LhsEval>(fluidState.massFraction(FluidSystem::oilPhaseIdx, FluidSystem::gasCompIdx));
    return FluidSystem::convertXoGToRs(XoG, regionIdx);
}
 
template <class FluidSystem, class FluidState, class LhsEval>
auto getRs_(typename std::enable_if<HasMember_Rs<FluidState>::value, const FluidState&>::type fluidState,
            unsigned)
    -> decltype(decay<LhsEval>(fluidState.Rs()))
{ return decay<LhsEval>(fluidState.Rs()); }
 
template <class FluidSystem, class FluidState, class LhsEval>
LhsEval getRv_(typename std::enable_if<!HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
               unsigned regionIdx)
{
    const auto& XgO =
        decay<LhsEval>(fluidState.massFraction(FluidSystem::gasPhaseIdx, FluidSystem::oilCompIdx));
    return FluidSystem::convertXgOToRv(XgO, regionIdx);
}
 
template <class FluidSystem, class FluidState, class LhsEval>
auto getRv_(typename std::enable_if<HasMember_Rv<FluidState>::value, const FluidState&>::type fluidState,
            unsigned)
    -> decltype(decay<LhsEval>(fluidState.Rv()))
{ return decay<LhsEval>(fluidState.Rv()); }
 
template <class FluidSystem, class FluidState, class LhsEval>
LhsEval getRvw_(typename std::enable_if<!HasMember_Rvw<FluidState>::value, const FluidState&>::type fluidState,
               unsigned regionIdx)
{
    const auto& XgW =
        decay<LhsEval>(fluidState.massFraction(FluidSystem::gasPhaseIdx, FluidSystem::waterCompIdx));
    return FluidSystem::convertXgWToRvw(XgW, regionIdx);
}
 
template <class FluidSystem, class FluidState, class LhsEval>
auto getRvw_(typename std::enable_if<HasMember_Rvw<FluidState>::value, const FluidState&>::type fluidState,
            unsigned)
    -> decltype(decay<LhsEval>(fluidState.Rvw()))
{ return decay<LhsEval>(fluidState.Rvw()); }
 
template <class FluidSystem, class FluidState, class LhsEval>
LhsEval getSaltConcentration_(typename std::enable_if<!HasMember_saltConcentration<FluidState>::value,
                              const FluidState&>::type,
                              unsigned)
{return 0.0;}
 
template <class FluidSystem, class FluidState, class LhsEval>
auto getSaltConcentration_(typename std::enable_if<HasMember_saltConcentration<FluidState>::value, const FluidState&>::type fluidState,
            unsigned)
    -> decltype(decay<LhsEval>(fluidState.saltConcentration()))
{ return decay<LhsEval>(fluidState.saltConcentration()); }
 
template <class FluidSystem, class FluidState, class LhsEval>
LhsEval getSaltSaturation_(typename std::enable_if<!HasMember_saltSaturation<FluidState>::value,
                              const FluidState&>::type,
                              unsigned)
{return 0.0;}
 
template <class FluidSystem, class FluidState, class LhsEval>
auto getSaltSaturation_(typename std::enable_if<HasMember_saltSaturation<FluidState>::value, const FluidState&>::type fluidState,
            unsigned)
    -> decltype(decay<LhsEval>(fluidState.saltSaturation()))
{ return decay<LhsEval>(fluidState.saltSaturation()); }
 
}
 
template <class Scalar, class IndexTraits = BlackOilDefaultIndexTraits>
class BlackOilFluidSystem : public BaseFluidSystem<Scalar, BlackOilFluidSystem<Scalar, IndexTraits> >
{
    using ThisType = BlackOilFluidSystem;
 
public:
    using GasPvt = GasPvtMultiplexer<Scalar>;
    using OilPvt = OilPvtMultiplexer<Scalar>;
    using WaterPvt = WaterPvtMultiplexer<Scalar>;
 
    template <class EvaluationT>
    struct ParameterCache : public NullParameterCache<EvaluationT>
    {
        using Evaluation = EvaluationT;
 
    public:
        ParameterCache(Scalar maxOilSat = 1.0, unsigned regionIdx=0)
        {
            maxOilSat_ = maxOilSat;
            regionIdx_ = regionIdx;
        }
 
        template <class OtherCache>
        void assignPersistentData(const OtherCache& other)
        {
            regionIdx_ = other.regionIndex();
            maxOilSat_ = other.maxOilSat();
        }
 
        unsigned regionIndex() const
        { return regionIdx_; }
 
        void setRegionIndex(unsigned val)
        { regionIdx_ = val; }
 
        const Evaluation& maxOilSat() const
        { return maxOilSat_; }
 
        void setMaxOilSat(const Evaluation& val)
        { maxOilSat_ = val; }
 
    private:
        Evaluation maxOilSat_;
        unsigned regionIdx_;
    };
 
    /****************************************
     * Initialization
     ****************************************/
#if HAVE_ECL_INPUT
    static void initFromState(const EclipseState& eclState, const Schedule& schedule)
    {
        size_t numRegions = eclState.runspec().tabdims().getNumPVTTables();
        initBegin(numRegions);
 
        numActivePhases_ = 0;
        std::fill_n(&phaseIsActive_[0], numPhases, false);
 
 
        if (eclState.runspec().phases().active(Phase::OIL)) {
            phaseIsActive_[oilPhaseIdx] = true;
            ++ numActivePhases_;
        }
 
        if (eclState.runspec().phases().active(Phase::GAS)) {
            phaseIsActive_[gasPhaseIdx] = true;
            ++ numActivePhases_;
        }
 
        if (eclState.runspec().phases().active(Phase::WATER)) {
            phaseIsActive_[waterPhaseIdx] = true;
            ++ numActivePhases_;
        }
 
        // set the surface conditions using the STCOND keyword
        surfaceTemperature = eclState.getTableManager().stCond().temperature;
        surfacePressure = eclState.getTableManager().stCond().pressure;
 
        // The reservoir temperature does not really belong into the table manager. TODO:
        // change this in opm-parser
        setReservoirTemperature(eclState.getTableManager().rtemp());
 
        // this fluidsystem only supports two or three phases
        assert(numActivePhases_ >= 1 && numActivePhases_ <= 3);
 
        setEnableDissolvedGas(eclState.getSimulationConfig().hasDISGAS());
        setEnableVaporizedOil(eclState.getSimulationConfig().hasVAPOIL());
        setEnableVaporizedWater(eclState.getSimulationConfig().hasVAPWAT());
 
        if (phaseIsActive(gasPhaseIdx)) {
            gasPvt_ = std::make_shared<GasPvt>();
            gasPvt_->initFromState(eclState, schedule);
        }
 
        if (phaseIsActive(oilPhaseIdx)) {
            oilPvt_ = std::make_shared<OilPvt>();
            oilPvt_->initFromState(eclState, schedule);
        }
 
        if (phaseIsActive(waterPhaseIdx)) {
            waterPvt_ = std::make_shared<WaterPvt>();
            waterPvt_->initFromState(eclState, schedule);
        }
 
        // set the reference densities of all PVT regions
        for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
            setReferenceDensities(phaseIsActive(oilPhaseIdx)? oilPvt_->oilReferenceDensity(regionIdx):700.,
                                  phaseIsActive(waterPhaseIdx)? waterPvt_->waterReferenceDensity(regionIdx):1000.,
                                  phaseIsActive(gasPhaseIdx)? gasPvt_->gasReferenceDensity(regionIdx):2.,
                                  regionIdx);
        }
 
        // set default molarMass and mappings
        initEnd();
 
        // use molarMass of CO2 and Brine as default
        // when we are using the the CO2STORE option
        // NB the oil component is used internally for
        // brine
        if (eclState.runspec().co2Storage()) {
            for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
                molarMass_[regionIdx][oilCompIdx] = BrineCo2Pvt<Scalar>::Brine::molarMass();
                molarMass_[regionIdx][gasCompIdx] = BrineCo2Pvt<Scalar>::CO2::molarMass();
            }
        }
 
        setEnableDiffusion(eclState.getSimulationConfig().isDiffusive());
        if(enableDiffusion()) {
            const auto& diffCoeffTables = eclState.getTableManager().getDiffusionCoefficientTable();
            if(!diffCoeffTables.empty()) {
                // if diffusion coefficient table is empty we relay on the PVT model to
                // to give us the coefficients.
                diffusionCoefficients_.resize(numRegions,{0,0,0,0,0,0,0,0,0});
                assert(diffCoeffTables.size() == numRegions);
                for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
                    const auto& diffCoeffTable = diffCoeffTables[regionIdx];
                    molarMass_[regionIdx][oilCompIdx] = diffCoeffTable.oil_mw;
                    molarMass_[regionIdx][gasCompIdx] = diffCoeffTable.gas_mw;
                    setDiffusionCoefficient(diffCoeffTable.gas_in_gas, gasCompIdx, gasPhaseIdx, regionIdx);
                    setDiffusionCoefficient(diffCoeffTable.oil_in_gas, oilCompIdx, gasPhaseIdx, regionIdx);
                    setDiffusionCoefficient(diffCoeffTable.gas_in_oil, gasCompIdx, oilPhaseIdx, regionIdx);
                    setDiffusionCoefficient(diffCoeffTable.oil_in_oil, oilCompIdx, oilPhaseIdx, regionIdx);
                    if(diffCoeffTable.gas_in_oil_cross_phase > 0 || diffCoeffTable.oil_in_oil_cross_phase > 0) {
                        throw std::runtime_error("Cross phase diffusion is set in the deck, but not implemented in Flow. "
                                                 "Please default DIFFC item 7 and item 8 or set it to zero.");
                    }
                }
            }
        }
    }
#endif // HAVE_ECL_INPUT
 
    static void initBegin(size_t numPvtRegions)
    {
        isInitialized_ = false;
 
        enableDissolvedGas_ = true;
        enableVaporizedOil_ = false;
        enableVaporizedWater_ = false;
        enableDiffusion_ = false;
 
        oilPvt_ = nullptr;
        gasPvt_ = nullptr;
        waterPvt_ = nullptr;
 
        surfaceTemperature = 273.15 + 15.56; // [K]
        surfacePressure = 1.01325e5; // [Pa]
        setReservoirTemperature(surfaceTemperature);
 
        numActivePhases_ = numPhases;
        std::fill_n(&phaseIsActive_[0], numPhases, true);
 
        resizeArrays_(numPvtRegions);
    }
 
    static void setEnableDissolvedGas(bool yesno)
    { enableDissolvedGas_ = yesno; }
 
    static void setEnableVaporizedOil(bool yesno)
    { enableVaporizedOil_ = yesno; }
 
    static void setEnableVaporizedWater(bool yesno)
    { enableVaporizedWater_ = yesno; }
 
    static void setEnableDiffusion(bool yesno)
    { enableDiffusion_ = yesno; }
 
 
    static void setGasPvt(std::shared_ptr<GasPvt> pvtObj)
    { gasPvt_ = pvtObj; }
 
    static void setOilPvt(std::shared_ptr<OilPvt> pvtObj)
    { oilPvt_ = pvtObj; }
 
    static void setWaterPvt(std::shared_ptr<WaterPvt> pvtObj)
    { waterPvt_ = pvtObj; }
 
    static void setReferenceDensities(Scalar rhoOil,
                                      Scalar rhoWater,
                                      Scalar rhoGas,
                                      unsigned regionIdx)
    {
        referenceDensity_[regionIdx][oilPhaseIdx] = rhoOil;
        referenceDensity_[regionIdx][waterPhaseIdx] = rhoWater;
        referenceDensity_[regionIdx][gasPhaseIdx] = rhoGas;
    }
 
 
    static void initEnd()
    {
        // calculate the final 2D functions which are used for interpolation.
        size_t numRegions = molarMass_.size();
        for (unsigned regionIdx = 0; regionIdx < numRegions; ++ regionIdx) {
            // calculate molar masses
 
            // water is simple: 18 g/mol
            molarMass_[regionIdx][waterCompIdx] = 18e-3;
 
            if (phaseIsActive(gasPhaseIdx)) {
                // for gas, we take the density at standard conditions and assume it to be ideal
                Scalar p = surfacePressure;
                Scalar T = surfaceTemperature;
                Scalar rho_g = referenceDensity_[/*regionIdx=*/0][gasPhaseIdx];
                molarMass_[regionIdx][gasCompIdx] = Constants<Scalar>::R*T*rho_g / p;
            }
            else
                // hydrogen gas. we just set this do avoid NaNs later
                molarMass_[regionIdx][gasCompIdx] = 2e-3;
 
            // finally, for oil phase, we take the molar mass from the spe9 paper
            molarMass_[regionIdx][oilCompIdx] = 175e-3; // kg/mol
        }
 
 
        int activePhaseIdx = 0;
        for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
            if(phaseIsActive(phaseIdx)){
                canonicalToActivePhaseIdx_[phaseIdx] = activePhaseIdx;
                activeToCanonicalPhaseIdx_[activePhaseIdx] = phaseIdx;
                activePhaseIdx++;
            }
        }
        isInitialized_ = true;
    }
 
    static bool isInitialized()
    { return isInitialized_; }
 
    /****************************************
     * Generic phase properties
     ****************************************/
 
    static constexpr unsigned numPhases = 3;
 
    static constexpr unsigned waterPhaseIdx = IndexTraits::waterPhaseIdx;
    static constexpr unsigned oilPhaseIdx = IndexTraits::oilPhaseIdx;
    static constexpr unsigned gasPhaseIdx = IndexTraits::gasPhaseIdx;
 
    static Scalar surfacePressure;
 
    static Scalar surfaceTemperature;
 
    static const char* phaseName(unsigned phaseIdx)
    {
        switch (phaseIdx) {
        case waterPhaseIdx:
            return "water";
        case oilPhaseIdx:
            return "oil";
        case gasPhaseIdx:
            return "gas";
 
        default:
            throw std::logic_error("Phase index " + std::to_string(phaseIdx) + " is unknown");
        }
    }
 
    static bool isLiquid(unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);
        return phaseIdx != gasPhaseIdx;
    }
 
    /****************************************
     * Generic component related properties
     ****************************************/
 
    static constexpr unsigned numComponents = 3;
 
    static constexpr unsigned oilCompIdx = IndexTraits::oilCompIdx;
    static constexpr unsigned waterCompIdx = IndexTraits::waterCompIdx;
    static constexpr unsigned gasCompIdx = IndexTraits::gasCompIdx;
 
protected:
    static unsigned char numActivePhases_;
    static std::array<bool,numPhases> phaseIsActive_;
 
public:
    static unsigned numActivePhases()
    { return numActivePhases_; }
 
    static unsigned phaseIsActive(unsigned phaseIdx)
    {
        assert(phaseIdx < numPhases);
        return phaseIsActive_[phaseIdx];
    }
 
    static constexpr unsigned solventComponentIndex(unsigned phaseIdx)
    {
        switch (phaseIdx) {
        case waterPhaseIdx:
            return waterCompIdx;
        case oilPhaseIdx:
            return oilCompIdx;
        case gasPhaseIdx:
            return gasCompIdx;
 
        default:
            throw std::logic_error("Phase index " + std::to_string(phaseIdx) + " is unknown");
        }
    }
 
    static constexpr unsigned soluteComponentIndex(unsigned phaseIdx)
    {
        switch (phaseIdx) {
        case waterPhaseIdx:
            throw std::logic_error("The water phase does not have any solutes in the black oil model!");
        case oilPhaseIdx:
            return gasCompIdx;
        case gasPhaseIdx:
            return oilCompIdx;
 
        default:
            throw std::logic_error("Phase index " + std::to_string(phaseIdx) + " is unknown");
        }
    }
 
    static const char* componentName(unsigned compIdx)
    {
        switch (compIdx) {
        case waterCompIdx:
            return "Water";
        case oilCompIdx:
            return "Oil";
        case gasCompIdx:
            return "Gas";
 
        default:
            throw std::logic_error("Component index " + std::to_string(compIdx) + " is unknown");
        }
    }
 
    static Scalar molarMass(unsigned compIdx, unsigned regionIdx = 0)
    { return molarMass_[regionIdx][compIdx]; }
 
    static bool isIdealMixture(unsigned /*phaseIdx*/)
    {
        // fugacity coefficients are only pressure dependent -> we
        // have an ideal mixture
        return true;
    }
 
    static bool isCompressible(unsigned /*phaseIdx*/)
    { return true; /* all phases are compressible */ }
 
    static bool isIdealGas(unsigned /*phaseIdx*/)
    { return false; }
 
 
    /****************************************
     * Black-oil specific properties
     ****************************************/
    static size_t numRegions()
    { return molarMass_.size(); }
 
    static bool enableDissolvedGas()
    { return enableDissolvedGas_; }
 
    static bool enableVaporizedOil()
    { return enableVaporizedOil_; }
 
    static bool enableVaporizedWater()
    { return enableVaporizedWater_; }
 
    static bool enableDiffusion()
    { return enableDiffusion_; }
 
    static Scalar referenceDensity(unsigned phaseIdx, unsigned regionIdx)
    { return referenceDensity_[regionIdx][phaseIdx]; }
 
    /****************************************
     * thermodynamic quantities (generic version)
     ****************************************/
    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
    static LhsEval density(const FluidState& fluidState,
                           const ParameterCache<ParamCacheEval>& paramCache,
                           unsigned phaseIdx)
    { return density<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
    static LhsEval fugacityCoefficient(const FluidState& fluidState,
                                       const ParameterCache<ParamCacheEval>& paramCache,
                                       unsigned phaseIdx,
                                       unsigned compIdx)
    {
        return fugacityCoefficient<FluidState, LhsEval>(fluidState,
                                                        phaseIdx,
                                                        compIdx,
                                                        paramCache.regionIndex());
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
    static LhsEval viscosity(const FluidState& fluidState,
                             const ParameterCache<ParamCacheEval>& paramCache,
                             unsigned phaseIdx)
    { return viscosity<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
    static LhsEval enthalpy(const FluidState& fluidState,
                            const ParameterCache<ParamCacheEval>& paramCache,
                            unsigned phaseIdx)
    { return enthalpy<FluidState, LhsEval>(fluidState, phaseIdx, paramCache.regionIndex()); }
 
    /****************************************
     * thermodynamic quantities (black-oil specific version: Note that the PVT region
     * index is explicitly passed instead of a parameter cache object)
     ****************************************/
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval density(const FluidState& fluidState,
                           unsigned phaseIdx,
                           unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& saltConcentration = BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
 
        switch (phaseIdx) {
        case oilPhaseIdx: {
            if (enableDissolvedGas()) {
                // miscible oil
                const LhsEval& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
 
                return
                    bo*referenceDensity(oilPhaseIdx, regionIdx)
                    + Rs*bo*referenceDensity(gasPhaseIdx, regionIdx);
            }
 
            // immiscible oil
            const LhsEval Rs(0.0);
            const auto& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
 
            return referenceDensity(phaseIdx, regionIdx)*bo;
        }
 
        case gasPhaseIdx: {
             if (enableVaporizedOil() && enableVaporizedWater()) {
                // gas containing vaporized oil and vaporized water
                const LhsEval& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                const LhsEval& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
                return
                    bg*referenceDensity(gasPhaseIdx, regionIdx)
                    + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx)
                    + Rvw*bg*referenceDensity(waterPhaseIdx, regionIdx);
            }
            if (enableVaporizedOil()) {
                // miscible gas
                const LhsEval Rvw(0.0);
                const LhsEval& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
                return
                    bg*referenceDensity(gasPhaseIdx, regionIdx)
                    + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
            }
            if (enableVaporizedWater()) {
                // gas containing vaporized water
                const LhsEval Rv(0.0);
                const LhsEval& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
                return
                    bg*referenceDensity(gasPhaseIdx, regionIdx)
                    + Rvw*bg*referenceDensity(waterPhaseIdx, regionIdx);
            }
 
            // immiscible gas
            const LhsEval Rv(0.0);
            const LhsEval Rvw(0.0);
            const auto& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
            return bg*referenceDensity(phaseIdx, regionIdx);
        }
 
        case waterPhaseIdx:
            return
                referenceDensity(waterPhaseIdx, regionIdx)
                * waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
        }
 
        throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval saturatedDensity(const FluidState& fluidState,
                                    unsigned phaseIdx,
                                    unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = fluidState.pressure(phaseIdx);
        const auto& T = fluidState.temperature(phaseIdx);
 
        switch (phaseIdx) {
        case oilPhaseIdx: {
            if (enableDissolvedGas()) {
                // miscible oil
                const LhsEval& Rs = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, oilPhaseIdx, regionIdx);
                const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
 
                return
                    bo*referenceDensity(oilPhaseIdx, regionIdx)
                    + Rs*bo*referenceDensity(gasPhaseIdx, regionIdx);
            }
 
            // immiscible oil
            const LhsEval Rs(0.0);
            const LhsEval& bo = oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
            return referenceDensity(phaseIdx, regionIdx)*bo;
        }
 
        case gasPhaseIdx: {
            if (enableVaporizedOil() && enableVaporizedWater()) {
                // gas containing vaporized oil and vaporized water
                const LhsEval& Rv = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
                const LhsEval& Rvw = saturatedVaporizationFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
                const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
                return
                    bg*referenceDensity(gasPhaseIdx, regionIdx)
                    + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx) 
                    + Rvw*bg*referenceDensity(waterPhaseIdx, regionIdx) ;
            }
 
            if (enableVaporizedOil()) {
                // miscible gas
                const LhsEval Rvw(0.0);
                const LhsEval& Rv = saturatedDissolutionFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
                const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
                return
                    bg*referenceDensity(gasPhaseIdx, regionIdx)
                    + Rv*bg*referenceDensity(oilPhaseIdx, regionIdx);
            }
 
            if (enableVaporizedWater()) {
                // gas containing vaporized water
                const LhsEval Rv(0.0);
                const LhsEval& Rvw = saturatedVaporizationFactor<FluidState, LhsEval>(fluidState, gasPhaseIdx, regionIdx);
                const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
                return
                    bg*referenceDensity(gasPhaseIdx, regionIdx)
                    + Rvw*bg*referenceDensity(waterPhaseIdx, regionIdx);
            }
 
            // immiscible gas
            const LhsEval Rv(0.0);
            const LhsEval Rvw(0.0);
            const LhsEval& bg = gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
 
            return referenceDensity(phaseIdx, regionIdx)*bg;
 
        }
 
        case waterPhaseIdx:
            return
                referenceDensity(waterPhaseIdx, regionIdx)
                *inverseFormationVolumeFactor<FluidState, LhsEval>(fluidState, waterPhaseIdx, regionIdx);
        }
 
        throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval inverseFormationVolumeFactor(const FluidState& fluidState,
                                                unsigned phaseIdx,
                                                unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& saltConcentration = decay<LhsEval>(fluidState.saltConcentration());
 
        switch (phaseIdx) {
        case oilPhaseIdx: {
            if (enableDissolvedGas()) {
                const auto& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                if (fluidState.saturation(gasPhaseIdx) > 0.0
                    && Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p)))
                {
                    return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
                } else {
                    return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
                }
            }
 
            const LhsEval Rs(0.0);
            return oilPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rs);
        }
        case gasPhaseIdx: {
            if (enableVaporizedOil() && enableVaporizedWater()) {
                 const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                 const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                 if (fluidState.saturation(waterPhaseIdx) > 0.0
                    && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p))
                    && fluidState.saturation(oilPhaseIdx) > 0.0
                    && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 { 
                    return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
                 } else {
                     return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
                 }
            }
 
            if (enableVaporizedOil()) {
                const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                if (fluidState.saturation(oilPhaseIdx) > 0.0
                    && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                {
                    return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
                } else {
                    const LhsEval Rvw(0.0);
                    return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
                }
            }
 
            if (enableVaporizedWater()) { 
                const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                if (fluidState.saturation(waterPhaseIdx) > 0.0
                    && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                {
                    return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
                } else {
                    const LhsEval Rv(0.0);
                    return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
                }
            }
            
            const LhsEval Rv(0.0);
            const LhsEval Rvw(0.0);
            return gasPvt_->inverseFormationVolumeFactor(regionIdx, T, p, Rv, Rvw);
        }
        case waterPhaseIdx:
            return waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval saturatedInverseFormationVolumeFactor(const FluidState& fluidState,
                                                         unsigned phaseIdx,
                                                         unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& saltConcentration = decay<LhsEval>(fluidState.saltConcentration());
 
        switch (phaseIdx) {
        case oilPhaseIdx: return oilPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
        case gasPhaseIdx: return gasPvt_->saturatedInverseFormationVolumeFactor(regionIdx, T, p);
        case waterPhaseIdx: return waterPvt_->inverseFormationVolumeFactor(regionIdx, T, p, saltConcentration);
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval fugacityCoefficient(const FluidState& fluidState,
                                       unsigned phaseIdx,
                                       unsigned compIdx,
                                       unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(compIdx <= numComponents);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
 
        // for the fugacity coefficient of the oil component in the oil phase, we use
        // some pseudo-realistic value for the vapor pressure to ease physical
        // interpretation of the results
        const LhsEval phi_oO = 20e3/p;
 
        // for the gas component in the gas phase, assume it to be an ideal gas
        constexpr const Scalar phi_gG = 1.0;
 
        // for the fugacity coefficient of the water component in the water phase, we use
        // the same approach as for the oil component in the oil phase
        const LhsEval phi_wW = 30e3/p;
 
        switch (phaseIdx) {
        case gasPhaseIdx: // fugacity coefficients for all components in the gas phase
            switch (compIdx) {
            case gasCompIdx:
                return phi_gG;
 
            // for the oil component, we calculate the Rv value for saturated gas and Rs
            // for saturated oil, and compute the fugacity coefficients at the
            // equilibrium. for this, we assume that in equilibrium the fugacities of the
            // oil component is the same in both phases.
            case oilCompIdx: {
                if (!enableVaporizedOil())
                    // if there's no vaporized oil, the gas phase is assumed to be
                    // immiscible with the oil component
                    return phi_gG*1e6;
 
                const auto& R_vSat = gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
                const auto& X_gOSat = convertRvToXgO(R_vSat, regionIdx);
                const auto& x_gOSat = convertXgOToxgO(X_gOSat, regionIdx);
 
                const auto& R_sSat = oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
                const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
                const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
                const auto& x_oOSat = 1.0 - x_oGSat;
 
                const auto& p_o = decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
                const auto& p_g = decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
 
                return phi_oO*p_o*x_oOSat / (p_g*x_gOSat);
            }
 
            case waterCompIdx:
                // the water component is assumed to be never miscible with the gas phase
                return phi_gG*1e6;
 
            default:
                throw std::logic_error("Invalid component index "+std::to_string(compIdx));
            }
 
        case oilPhaseIdx: // fugacity coefficients for all components in the oil phase
            switch (compIdx) {
            case oilCompIdx:
                return phi_oO;
 
            // for the oil and water components, we proceed analogous to the gas and
            // water components in the gas phase
            case gasCompIdx: {
                if (!enableDissolvedGas())
                    // if there's no dissolved gas, the oil phase is assumed to be
                    // immiscible with the gas component
                    return phi_oO*1e6;
 
                const auto& R_vSat = gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
                const auto& X_gOSat = convertRvToXgO(R_vSat, regionIdx);
                const auto& x_gOSat = convertXgOToxgO(X_gOSat, regionIdx);
                const auto& x_gGSat = 1.0 - x_gOSat;
 
                const auto& R_sSat = oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
                const auto& X_oGSat = convertRsToXoG(R_sSat, regionIdx);
                const auto& x_oGSat = convertXoGToxoG(X_oGSat, regionIdx);
 
                const auto& p_o = decay<LhsEval>(fluidState.pressure(oilPhaseIdx));
                const auto& p_g = decay<LhsEval>(fluidState.pressure(gasPhaseIdx));
 
                return phi_gG*p_g*x_gGSat / (p_o*x_oGSat);
            }
 
            case waterCompIdx:
                return phi_oO*1e6;
 
            default:
                throw std::logic_error("Invalid component index "+std::to_string(compIdx));
            }
 
        case waterPhaseIdx: // fugacity coefficients for all components in the water phase
            // the water phase fugacity coefficients are pretty simple: because the water
            // phase is assumed to consist entirely from the water component, we just
            // need to make sure that the fugacity coefficients for the other components
            // are a few orders of magnitude larger than the one of the water
            // component. (i.e., the affinity of the gas and oil components to the water
            // phase is lower by a few orders of magnitude)
            switch (compIdx) {
            case waterCompIdx: return phi_wW;
            case oilCompIdx: return 1.1e6*phi_wW;
            case gasCompIdx: return 1e6*phi_wW;
            default:
                throw std::logic_error("Invalid component index "+std::to_string(compIdx));
            }
 
        default:
            throw std::logic_error("Invalid phase index "+std::to_string(phaseIdx));
        }
 
        throw std::logic_error("Unhandled phase or component index");
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval viscosity(const FluidState& fluidState,
                             unsigned phaseIdx,
                             unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const LhsEval& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const LhsEval& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const LhsEval& saltConcentration = BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
 
        switch (phaseIdx) {
        case oilPhaseIdx: {
            if (enableDissolvedGas()) {
                const auto& Rs = BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                if (fluidState.saturation(gasPhaseIdx) > 0.0
                    && Rs >= (1.0 - 1e-10)*oilPvt_->saturatedGasDissolutionFactor(regionIdx, scalarValue(T), scalarValue(p)))
                {
                    return oilPvt_->saturatedViscosity(regionIdx, T, p);
                } else {
                    return oilPvt_->viscosity(regionIdx, T, p, Rs);
                }
            }
 
            const LhsEval Rs(0.0);
            return oilPvt_->viscosity(regionIdx, T, p, Rs);
        }
 
        case gasPhaseIdx: {
             if (enableVaporizedOil() && enableVaporizedWater()) {
                 const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                 const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                 if (fluidState.saturation(waterPhaseIdx) > 0.0
                    && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p))
                    && fluidState.saturation(oilPhaseIdx) > 0.0
                    && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                 { 
                     return gasPvt_->saturatedViscosity(regionIdx, T, p);
                 } else {
                     return gasPvt_->viscosity(regionIdx, T, p, Rv, Rvw);
                 }
            }
            if (enableVaporizedOil()) {
                const auto& Rv = BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                if (fluidState.saturation(oilPhaseIdx) > 0.0
                    && Rv >= (1.0 - 1e-10)*gasPvt_->saturatedOilVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                {
                    return gasPvt_->saturatedViscosity(regionIdx, T, p);
                } else {
                    const LhsEval Rvw(0.0);
                    return gasPvt_->viscosity(regionIdx, T, p, Rv, Rvw);
                }
            }
            if (enableVaporizedWater()) {
                const auto& Rvw = BlackOil::template getRvw_<ThisType, FluidState, LhsEval>(fluidState, regionIdx);
                if (fluidState.saturation(waterPhaseIdx) > 0.0
                    && Rvw >= (1.0 - 1e-10)*gasPvt_->saturatedWaterVaporizationFactor(regionIdx, scalarValue(T), scalarValue(p)))
                {
                    return gasPvt_->saturatedViscosity(regionIdx, T, p); 
                } else {
                    const LhsEval Rv(0.0);
                    return gasPvt_->viscosity(regionIdx, T, p, Rv, Rvw);
                }
            }
 
            const LhsEval Rv(0.0);
            const LhsEval Rvw(0.0);
            return gasPvt_->viscosity(regionIdx, T, p, Rv, Rvw);
        }
 
        case waterPhaseIdx:
            // since water is always assumed to be immiscible in the black-oil model,
            // there is no "saturated water"
            return waterPvt_->viscosity(regionIdx, T, p, saltConcentration);
        }
 
        throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval enthalpy(const FluidState& fluidState,
                            unsigned phaseIdx,
                            unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
 
        switch (phaseIdx) {
        case oilPhaseIdx:
            return
                oilPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
                + p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
 
        case gasPhaseIdx:
            return
                 gasPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
                  + p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
 
        case waterPhaseIdx:
            return
                waterPvt_->internalEnergy(regionIdx, T, p, BlackOil::template getSaltConcentration_<ThisType, FluidState, LhsEval>(fluidState, regionIdx))
                + p/density<FluidState, LhsEval>(fluidState, phaseIdx, regionIdx);
 
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
 
        throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval saturatedVaporizationFactor(const FluidState& fluidState,
                                              unsigned phaseIdx,
                                              unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& saltConcentration = decay<LhsEval>(fluidState.saltConcentration());
 
        switch (phaseIdx) {
        case oilPhaseIdx: return 0.0;
        case gasPhaseIdx: return gasPvt_->saturatedWaterVaporizationFactor(regionIdx, T, p, saltConcentration);
        case waterPhaseIdx: return 0.0;
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval saturatedDissolutionFactor(const FluidState& fluidState,
                                              unsigned phaseIdx,
                                              unsigned regionIdx,
                                              const LhsEval& maxOilSaturation)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
        const auto& So = decay<LhsEval>(fluidState.saturation(oilPhaseIdx));
 
        switch (phaseIdx) {
        case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p, So, maxOilSaturation);
        case gasPhaseIdx: return gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p, So, maxOilSaturation);
        case waterPhaseIdx: return 0.0;
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval saturatedDissolutionFactor(const FluidState& fluidState,
                                              unsigned phaseIdx,
                                              unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
 
        switch (phaseIdx) {
        case oilPhaseIdx: return oilPvt_->saturatedGasDissolutionFactor(regionIdx, T, p);
        case gasPhaseIdx: return gasPvt_->saturatedOilVaporizationFactor(regionIdx, T, p);
        case waterPhaseIdx: return 0.0;
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval bubblePointPressure(const FluidState& fluidState,
                                       unsigned regionIdx)
    {
        return saturationPressure(fluidState, oilPhaseIdx, regionIdx);
    }
 
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval dewPointPressure(const FluidState& fluidState,
                                       unsigned regionIdx)
    {
        return saturationPressure(fluidState, gasPhaseIdx, regionIdx);
    }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar>
    static LhsEval saturationPressure(const FluidState& fluidState,
                                      unsigned phaseIdx,
                                      unsigned regionIdx)
    {
        assert(phaseIdx <= numPhases);
        assert(regionIdx <= numRegions());
 
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
 
        switch (phaseIdx) {
        case oilPhaseIdx: return oilPvt_->saturationPressure(regionIdx, T, BlackOil::template getRs_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
        case gasPhaseIdx: return gasPvt_->saturationPressure(regionIdx, T, BlackOil::template getRv_<ThisType, FluidState, LhsEval>(fluidState, regionIdx));
        case waterPhaseIdx: return 0.0;
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
    /****************************************
     * Auxiliary and convenience methods for the black-oil model
     ****************************************/
    template <class LhsEval>
    static LhsEval convertXoGToRs(const LhsEval& XoG, unsigned regionIdx)
    {
        Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
        Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
 
        return XoG/(1.0 - XoG)*(rho_oRef/rho_gRef);
    }
 
    template <class LhsEval>
    static LhsEval convertXgOToRv(const LhsEval& XgO, unsigned regionIdx)
    {
        Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
        Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
 
        return XgO/(1.0 - XgO)*(rho_gRef/rho_oRef);
    }
 
    template <class LhsEval>
    static LhsEval convertXgWToRvw(const LhsEval& XgW, unsigned regionIdx)
    {
        Scalar rho_wRef = referenceDensity_[regionIdx][waterPhaseIdx];
        Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
 
        return XgW/(1.0 - XgW)*(rho_gRef/rho_wRef);
    }
 
 
    template <class LhsEval>
    static LhsEval convertRsToXoG(const LhsEval& Rs, unsigned regionIdx)
    {
        Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
        Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
 
        const LhsEval& rho_oG = Rs*rho_gRef;
        return rho_oG/(rho_oRef + rho_oG);
    }
 
    template <class LhsEval>
    static LhsEval convertRvToXgO(const LhsEval& Rv, unsigned regionIdx)
    {
        Scalar rho_oRef = referenceDensity_[regionIdx][oilPhaseIdx];
        Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
 
        const LhsEval& rho_gO = Rv*rho_oRef;
        return rho_gO/(rho_gRef + rho_gO);
    }
 
    template <class LhsEval>
    static LhsEval convertRvwToXgW(const LhsEval& Rvw, unsigned regionIdx)
    {
        Scalar rho_wRef = referenceDensity_[regionIdx][waterPhaseIdx];
        Scalar rho_gRef = referenceDensity_[regionIdx][gasPhaseIdx];
 
        const LhsEval& rho_gW = Rvw*rho_wRef;
        return rho_gW/(rho_gRef + rho_gW);
    }
 
    template <class LhsEval>
    static LhsEval convertXoGToxoG(const LhsEval& XoG, unsigned regionIdx)
    {
        Scalar MO = molarMass_[regionIdx][oilCompIdx];
        Scalar MG = molarMass_[regionIdx][gasCompIdx];
 
        return XoG*MO / (MG*(1 - XoG) + XoG*MO);
    }
 
    template <class LhsEval>
    static LhsEval convertxoGToXoG(const LhsEval& xoG, unsigned regionIdx)
    {
        Scalar MO = molarMass_[regionIdx][oilCompIdx];
        Scalar MG = molarMass_[regionIdx][gasCompIdx];
 
        return xoG*MG / (xoG*(MG - MO) + MO);
    }
 
    template <class LhsEval>
    static LhsEval convertXgOToxgO(const LhsEval& XgO, unsigned regionIdx)
    {
        Scalar MO = molarMass_[regionIdx][oilCompIdx];
        Scalar MG = molarMass_[regionIdx][gasCompIdx];
 
        return XgO*MG / (MO*(1 - XgO) + XgO*MG);
    }
 
    template <class LhsEval>
    static LhsEval convertxgOToXgO(const LhsEval& xgO, unsigned regionIdx)
    {
        Scalar MO = molarMass_[regionIdx][oilCompIdx];
        Scalar MG = molarMass_[regionIdx][gasCompIdx];
 
        return xgO*MO / (xgO*(MO - MG) + MG);
    }
 
    static const GasPvt& gasPvt()
    { return *gasPvt_; }
 
    static const OilPvt& oilPvt()
    { return *oilPvt_; }
 
    static const WaterPvt& waterPvt()
    { return *waterPvt_; }
 
    static Scalar reservoirTemperature(unsigned = 0)
    { return reservoirTemperature_; }
 
    static void setReservoirTemperature(Scalar value)
    { reservoirTemperature_ = value; }
 
    static short activeToCanonicalPhaseIdx(unsigned activePhaseIdx) {
        assert(activePhaseIdx<numActivePhases());
        return activeToCanonicalPhaseIdx_[activePhaseIdx];
    }
 
    static short canonicalToActivePhaseIdx(unsigned phaseIdx) {
        assert(phaseIdx<numPhases);
        assert(phaseIsActive(phaseIdx));
        return canonicalToActivePhaseIdx_[phaseIdx];
    }
 
    static Scalar diffusionCoefficient(unsigned compIdx, unsigned phaseIdx, unsigned regionIdx = 0)
    { return diffusionCoefficients_[regionIdx][numPhases*compIdx + phaseIdx]; }
 
    static void setDiffusionCoefficient(Scalar coefficient, unsigned compIdx, unsigned phaseIdx, unsigned regionIdx = 0)
    { diffusionCoefficients_[regionIdx][numPhases*compIdx + phaseIdx] = coefficient ; }
 
    template <class FluidState, class LhsEval = typename FluidState::Scalar, class ParamCacheEval = LhsEval>
    static LhsEval diffusionCoefficient(const FluidState& fluidState,
                                        const ParameterCache<ParamCacheEval>& paramCache,
                                        unsigned phaseIdx,
                                        unsigned compIdx)
    {
        // diffusion is disabled by the user
        if(!enableDiffusion())
            return 0.0;
 
        // diffusion coefficients are set, and we use them
        if(!diffusionCoefficients_.empty()) {
            return diffusionCoefficient(compIdx, phaseIdx, paramCache.regionIndex());
        }
 
        const auto& p = decay<LhsEval>(fluidState.pressure(phaseIdx));
        const auto& T = decay<LhsEval>(fluidState.temperature(phaseIdx));
 
        switch (phaseIdx) {
        case oilPhaseIdx: return oilPvt().diffusionCoefficient(T, p, compIdx);
        case gasPhaseIdx: return gasPvt().diffusionCoefficient(T, p, compIdx);
        case waterPhaseIdx: return 0.0;
        default: throw std::logic_error("Unhandled phase index "+std::to_string(phaseIdx));
        }
    }
 
private:
    static void resizeArrays_(size_t numRegions)
    {
        molarMass_.resize(numRegions);
        referenceDensity_.resize(numRegions);
    }
 
    static Scalar reservoirTemperature_;
 
    static std::shared_ptr<GasPvt> gasPvt_;
    static std::shared_ptr<OilPvt> oilPvt_;
    static std::shared_ptr<WaterPvt> waterPvt_;
 
    static bool enableDissolvedGas_;
    static bool enableVaporizedOil_;
    static bool enableVaporizedWater_;
    static bool enableDiffusion_;
 
    // HACK for GCC 4.4: the array size has to be specified using the literal value '3'
    // here, because GCC 4.4 seems to be unable to determine the number of phases from
    // the BlackOil fluid system in the attribute declaration below...
    static std::vector<std::array<Scalar, /*numPhases=*/3> > referenceDensity_;
    static std::vector<std::array<Scalar, /*numComponents=*/3> > molarMass_;
    static std::vector<std::array<Scalar, /*numComponents=*/3 * /*numPhases=*/3> > diffusionCoefficients_;
 
    static std::array<short, numPhases> activeToCanonicalPhaseIdx_;
    static std::array<short, numPhases> canonicalToActivePhaseIdx_;
 
    static bool isInitialized_;
};
 
template <class Scalar, class IndexTraits>
unsigned char BlackOilFluidSystem<Scalar, IndexTraits>::numActivePhases_;
 
template <class Scalar, class IndexTraits>
std::array<bool, BlackOilFluidSystem<Scalar, IndexTraits>::numPhases> BlackOilFluidSystem<Scalar, IndexTraits>::phaseIsActive_;
 
template <class Scalar, class IndexTraits>
std::array<short, BlackOilFluidSystem<Scalar, IndexTraits>::numPhases> BlackOilFluidSystem<Scalar, IndexTraits>::activeToCanonicalPhaseIdx_;
 
template <class Scalar, class IndexTraits>
std::array<short, BlackOilFluidSystem<Scalar, IndexTraits>::numPhases> BlackOilFluidSystem<Scalar, IndexTraits>::canonicalToActivePhaseIdx_;
 
template <class Scalar, class IndexTraits>
Scalar
BlackOilFluidSystem<Scalar, IndexTraits>::surfaceTemperature; // [K]
 
template <class Scalar, class IndexTraits>
Scalar
BlackOilFluidSystem<Scalar, IndexTraits>::surfacePressure; // [Pa]
 
template <class Scalar, class IndexTraits>
Scalar
BlackOilFluidSystem<Scalar, IndexTraits>::reservoirTemperature_;
 
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::enableDissolvedGas_;
 
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::enableVaporizedOil_;
 
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::enableVaporizedWater_;
 
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::enableDiffusion_;
 
template <class Scalar, class IndexTraits>
std::shared_ptr<OilPvtMultiplexer<Scalar> >
BlackOilFluidSystem<Scalar, IndexTraits>::oilPvt_;
 
template <class Scalar, class IndexTraits>
std::shared_ptr<GasPvtMultiplexer<Scalar> >
BlackOilFluidSystem<Scalar, IndexTraits>::gasPvt_;
 
template <class Scalar, class IndexTraits>
std::shared_ptr<WaterPvtMultiplexer<Scalar> >
BlackOilFluidSystem<Scalar, IndexTraits>::waterPvt_;
 
template <class Scalar, class IndexTraits>
std::vector<std::array<Scalar, 3> >
BlackOilFluidSystem<Scalar, IndexTraits>::referenceDensity_;
 
template <class Scalar, class IndexTraits>
std::vector<std::array<Scalar, 3> >
BlackOilFluidSystem<Scalar, IndexTraits>::molarMass_;
 
template <class Scalar, class IndexTraits>
std::vector<std::array<Scalar, 9> >
BlackOilFluidSystem<Scalar, IndexTraits>::diffusionCoefficients_;
 
template <class Scalar, class IndexTraits>
bool BlackOilFluidSystem<Scalar, IndexTraits>::isInitialized_ = false;
 
} // namespace Opm
 
#endif