BlackOilFluidState.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_STATE_HH
#define OPM_BLACK_OIL_FLUID_STATE_HH
 
#include <type_traits>
 
#include <opm/common/utility/gpuDecorators.hpp>
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/common/HasMemberGeneratorMacros.hpp>
 
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/ConditionalStorage.hpp>
 
namespace Opm {
 
OPM_GENERATE_HAS_MEMBER(pvtRegionIndex, ) // Creates 'HasMember_pvtRegionIndex<T>'.
 
template <class FluidState>
OPM_HOST_DEVICE unsigned getPvtRegionIndex_(typename std::enable_if<HasMember_pvtRegionIndex<FluidState>::value,
                                                    const FluidState&>::type fluidState)
{ return fluidState.pvtRegionIndex(); }
 
template <class FluidState>
OPM_HOST_DEVICE unsigned getPvtRegionIndex_(typename std::enable_if<!HasMember_pvtRegionIndex<FluidState>::value,
                                                    const FluidState&>::type)
{ return 0; }
 
OPM_GENERATE_HAS_MEMBER(invB, /*phaseIdx=*/0) // Creates 'HasMember_invB<T>'.
 
template <class FluidSystem, class FluidState, class LhsEval>
OPM_HOST_DEVICE
auto getInvB_(typename std::enable_if<HasMember_invB<FluidState>::value,
                                      const FluidState&>::type fluidState,
              unsigned phaseIdx,
              unsigned,
              [[maybe_unused]] const FluidSystem& fluidSystem = FluidSystem{})
    -> decltype(decay<LhsEval>(fluidState.invB(phaseIdx)))
{ return decay<LhsEval>(fluidState.invB(phaseIdx)); }
 
template <class FluidSystem, class FluidState, class LhsEval>
OPM_HOST_DEVICE
LhsEval getInvB_(typename std::enable_if<!HasMember_invB<FluidState>::value,
                                          const FluidState&>::type fluidState,
                 unsigned phaseIdx,
                 unsigned pvtRegionIdx,
                 const FluidSystem& fluidSystem = FluidSystem{})
{
    const auto& rho = fluidState.density(phaseIdx);
    const auto& Xsolvent =
        fluidState.massFraction(phaseIdx, fluidSystem.solventComponentIndex(phaseIdx));
 
    return
        decay<LhsEval>(rho)
        *decay<LhsEval>(Xsolvent)
        /fluidSystem.referenceDensity(phaseIdx, pvtRegionIdx);
}
 
OPM_GENERATE_HAS_MEMBER(saltConcentration, ) // Creates 'HasMember_saltConcentration<T>'.
 
template <class FluidState>
OPM_HOST_DEVICE auto getSaltConcentration_(typename std::enable_if<HasMember_saltConcentration<FluidState>::value,
                                                    const FluidState&>::type fluidState)
{ return fluidState.saltConcentration(); }
 
template <class FluidState>
OPM_HOST_DEVICE auto getSaltConcentration_(typename std::enable_if<!HasMember_saltConcentration<FluidState>::value,
                                                    const FluidState&>::type)
{ return 0.0; }
 
OPM_GENERATE_HAS_MEMBER(saltSaturation, ) // Creates 'HasMember_saltSaturation<T>'.
 
template <class FluidState>
OPM_HOST_DEVICE auto getSaltSaturation_(typename std::enable_if<HasMember_saltSaturation<FluidState>::value,
                                                    const FluidState&>::type fluidState)
{ return fluidState.saltSaturation(); }
 
 
template <class FluidState>
OPM_HOST_DEVICE auto getSaltSaturation_(typename std::enable_if<!HasMember_saltSaturation<FluidState>::value,
                                                    const FluidState&>::type)
{ return 0.0; }

template <class ScalarT,
          class FluidSystemT,
          bool enableTemperature = false,
          bool enableEnergy = false,
          bool enableDissolution = true,
          bool enableVapwat = false,
          bool enableBrine = false,
          bool enableSaltPrecipitation = false,
          bool enableDissolutionInWater = false,
          unsigned numStoragePhases = FluidSystemT::numPhases>
class BlackOilFluidState
{
public:
    using FluidSystem = FluidSystemT;
    using Scalar = ScalarT;
 
    static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
    static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
    static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
 
    static constexpr int waterCompIdx = FluidSystem::waterCompIdx;
    static constexpr int gasCompIdx = FluidSystem::gasCompIdx;
    static constexpr int oilCompIdx = FluidSystem::oilCompIdx;
 
    static constexpr int numPhases = FluidSystem::numPhases;
    static constexpr int numComponents = FluidSystem::numComponents;
 
    static constexpr bool fluidSystemIsStatic = std::is_empty_v<FluidSystem>;

    explicit OPM_HOST_DEVICE BlackOilFluidState(const FluidSystem& fluidSystem)
    {
        if constexpr (fluidSystemIsStatic) {
            fluidSystemPtr_ = &fluidSystem;
        }
    }

    OPM_HOST_DEVICE BlackOilFluidState()
    {
        static_assert(fluidSystemIsStatic);
    }

    void checkDefined() const
    {
#ifndef NDEBUG
        Valgrind::CheckDefined(pvtRegionIdx_);
 
        for (unsigned storagePhaseIdx = 0; storagePhaseIdx < numStoragePhases; ++ storagePhaseIdx) {
            Valgrind::CheckDefined(saturation_[storagePhaseIdx]);
            Valgrind::CheckDefined(pressure_[storagePhaseIdx]);
            Valgrind::CheckDefined(density_[storagePhaseIdx]);
            Valgrind::CheckDefined(invB_[storagePhaseIdx]);
 
            if constexpr (enableEnergy)
                Valgrind::CheckDefined((*enthalpy_)[storagePhaseIdx]);
        }
 
        if constexpr (enableDissolution) {
            Valgrind::CheckDefined(*Rs_);
            Valgrind::CheckDefined(*Rv_);
        }
 
        if constexpr (enableVapwat) {
            Valgrind::CheckDefined(*Rvw_);
        }
 
        if constexpr (enableDissolutionInWater) {
            Valgrind::CheckDefined(*Rsw_);
        }
 
        if constexpr (enableBrine) {
            Valgrind::CheckDefined(*saltConcentration_);
        }
 
        if constexpr (enableSaltPrecipitation) {
            Valgrind::CheckDefined(*saltSaturation_);
        }
 
        if constexpr (enableTemperature || enableEnergy)
            Valgrind::CheckDefined(*temperature_);
#endif // NDEBUG
    }

    template <class FluidState>
    OPM_HOST_DEVICE void assign(const FluidState& fs)
    {
        if constexpr (enableTemperature || enableEnergy)
            setTemperature(fs.temperature(/*phaseIdx=*/0));
 
        unsigned pvtRegionIdx = getPvtRegionIndex_<FluidState>(fs);
        setPvtRegionIndex(pvtRegionIdx);
 
        if constexpr (enableDissolution) {
            setRs(BlackOil::getRs_<FluidSystem, FluidState, Scalar>(fs, pvtRegionIdx));
            setRv(BlackOil::getRv_<FluidSystem, FluidState, Scalar>(fs, pvtRegionIdx));
        }
        if constexpr (enableVapwat) {
            setRvw(BlackOil::getRvw_<FluidSystem, FluidState, Scalar>(fs, pvtRegionIdx));
        }
        if constexpr (enableDissolutionInWater) {
            setRsw(BlackOil::getRsw_<FluidSystem, FluidState, Scalar>(fs, pvtRegionIdx));
        }
        if constexpr (enableBrine){
            setSaltConcentration(BlackOil::getSaltConcentration_<FluidState, Scalar>(fs, pvtRegionIdx));
        }
        if constexpr (enableSaltPrecipitation){
            setSaltSaturation(BlackOil::getSaltSaturation_<FluidSystem, FluidState, Scalar>(fs, pvtRegionIdx));
        }
        for (unsigned storagePhaseIdx = 0; storagePhaseIdx < numStoragePhases; ++storagePhaseIdx) {
            unsigned phaseIdx = storageToCanonicalPhaseIndex_(storagePhaseIdx, fluidSystem());
            setSaturation(phaseIdx, fs.saturation(phaseIdx));
            setPressure(phaseIdx, fs.pressure(phaseIdx));
            setDensity(phaseIdx, fs.density(phaseIdx));
 
            if constexpr (enableEnergy)
                setEnthalpy(phaseIdx, fs.enthalpy(phaseIdx));
 
            setInvB(phaseIdx, getInvB_<FluidSystem, FluidState, Scalar>(fs, phaseIdx, pvtRegionIdx, fluidSystem()));
        }
    }

    OPM_HOST_DEVICE void setPvtRegionIndex(unsigned newPvtRegionIdx)
    { pvtRegionIdx_ = static_cast<unsigned short>(newPvtRegionIdx); }

    OPM_HOST_DEVICE void setPressure(unsigned phaseIdx, const Scalar& p)
    { pressure_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())] = p; }

    OPM_HOST_DEVICE void setSaturation(unsigned phaseIdx, const Scalar& S)
    { saturation_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())] = S; }

    OPM_HOST_DEVICE void setTotalSaturation(const Scalar& value)
    {
        totalSaturation_ = value;
    }

    OPM_HOST_DEVICE void setTemperature(const Scalar& value)
    {
        assert(enableTemperature || enableEnergy);
 
        (*temperature_) = value;
    }

    OPM_HOST_DEVICE void setEnthalpy(unsigned phaseIdx, const Scalar& value)
    {
        assert(enableTemperature || enableEnergy);
 
        (*enthalpy_)[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())] = value;
    }

    OPM_HOST_DEVICE void setInvB(unsigned phaseIdx, const Scalar& b)
    { invB_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())] = b; }

    OPM_HOST_DEVICE void setDensity(unsigned phaseIdx, const Scalar& rho)
    { density_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())] = rho; }

    OPM_HOST_DEVICE void setRs(const Scalar& newRs)
    { *Rs_ = newRs; }

    OPM_HOST_DEVICE void setRv(const Scalar& newRv)
    { *Rv_ = newRv; }

    OPM_HOST_DEVICE void setRvw(const Scalar& newRvw)
    { *Rvw_ = newRvw; }

    OPM_HOST_DEVICE void setRsw(const Scalar& newRsw)
    { *Rsw_ = newRsw; }

    OPM_HOST_DEVICE void setSaltConcentration(const Scalar& newSaltConcentration)
    { *saltConcentration_ = newSaltConcentration; }

    OPM_HOST_DEVICE void setSaltSaturation(const Scalar& newSaltSaturation)
    { *saltSaturation_ = newSaltSaturation; }

    OPM_HOST_DEVICE const Scalar& pressure(unsigned phaseIdx) const
    { return pressure_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())]; }

    OPM_HOST_DEVICE const Scalar& saturation(unsigned phaseIdx) const
    { return saturation_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())]; }

    OPM_HOST_DEVICE const Scalar& totalSaturation() const
    {
        return totalSaturation_;
    }

    OPM_HOST_DEVICE Scalar temperature(unsigned) const
    {
        if constexpr (enableTemperature || enableEnergy) {
            return *temperature_;
        } else {
            return fluidSystem().reservoirTemperature(pvtRegionIdx_);
        }
    }

    OPM_HOST_DEVICE const Scalar& invB(unsigned phaseIdx) const
    { return invB_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())]; }

    OPM_HOST_DEVICE Scalar Rs() const
    {
        if constexpr (enableDissolution) {
            return *Rs_;
        } else {
            return Scalar(0.0);
        }
    }

    OPM_HOST_DEVICE Scalar Rv() const
    {
        if constexpr (!enableDissolution) {
            return Scalar(0.0);
        } else {
            return *Rv_;
        }
    }

    OPM_HOST_DEVICE Scalar Rvw() const
    {
        if constexpr (enableVapwat) {
            return *Rvw_;
        } else {
            return Scalar(0.0);
        }
    }

    OPM_HOST_DEVICE Scalar Rsw() const
    {
        if constexpr (enableDissolutionInWater) {
            return *Rsw_;
        } else {
            return Scalar(0.0);
        }
    }

    OPM_HOST_DEVICE Scalar saltConcentration() const
    {
        if constexpr (enableBrine) {
            return *saltConcentration_;
        } else {
            return Scalar(0.0);
        }
    }

    OPM_HOST_DEVICE Scalar saltSaturation() const
    {
        if constexpr (enableSaltPrecipitation) {
            return *saltSaturation_;
        } else {
            return Scalar(0.0);
        }
    }

    OPM_HOST_DEVICE unsigned short pvtRegionIndex() const
    { return pvtRegionIdx_; }

    OPM_HOST_DEVICE Scalar density(unsigned phaseIdx) const
    { return density_[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())]; }

    OPM_HOST_DEVICE const Scalar& enthalpy(unsigned phaseIdx) const
    { return (*enthalpy_)[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())]; }

    OPM_HOST_DEVICE Scalar internalEnergy(unsigned phaseIdx) const
    {   auto energy = (*enthalpy_)[canonicalToStoragePhaseIndex_(phaseIdx, fluidSystem())];
        if(!fluidSystem().enthalpyEqualEnergy()){
            energy -= pressure(phaseIdx)/density(phaseIdx);
        }
        return energy;
    }

    // slow methods

    OPM_HOST_DEVICE Scalar molarDensity(unsigned phaseIdx) const
    {
        const auto& rho = density(phaseIdx);
 
        if (phaseIdx == waterPhaseIdx)
            return rho/fluidSystem().molarMass(waterCompIdx, pvtRegionIdx_);
 
        return
            rho*(moleFraction(phaseIdx, gasCompIdx)/fluidSystem().molarMass(gasCompIdx, pvtRegionIdx_)
                 + moleFraction(phaseIdx, oilCompIdx)/fluidSystem().molarMass(oilCompIdx, pvtRegionIdx_));
 
    }

    OPM_HOST_DEVICE Scalar molarVolume(unsigned phaseIdx) const
    { return 1.0/molarDensity(phaseIdx); }

    OPM_HOST_DEVICE Scalar viscosity(unsigned phaseIdx) const
    { return fluidSystem().viscosity(*this, phaseIdx, pvtRegionIdx_); }

    OPM_HOST_DEVICE Scalar massFraction(unsigned phaseIdx, unsigned compIdx) const
    {
        switch (phaseIdx) {
        case waterPhaseIdx:
            if (compIdx == waterCompIdx)
                return 1.0;
            return 0.0;
 
        case oilPhaseIdx:
            if (compIdx == waterCompIdx)
                return 0.0;
            else if (compIdx == oilCompIdx)
                return 1.0 - fluidSystem().convertRsToXoG(Rs(), pvtRegionIdx_);
            else {
                assert(compIdx == gasCompIdx);
                return fluidSystem().convertRsToXoG(Rs(), pvtRegionIdx_);
            }
            break;
 
        case gasPhaseIdx:
            if (compIdx == waterCompIdx)
                return 0.0;
            else if (compIdx == oilCompIdx)
                return fluidSystem().convertRvToXgO(Rv(), pvtRegionIdx_);
            else {
                assert(compIdx == gasCompIdx);
                return 1.0 - fluidSystem().convertRvToXgO(Rv(), pvtRegionIdx_);
            }
            break;
        }
 
        OPM_THROW(std::logic_error, "Invalid phase or component index!");
    }

    OPM_HOST_DEVICE Scalar moleFraction(unsigned phaseIdx, unsigned compIdx) const
    {
        switch (phaseIdx) {
        case waterPhaseIdx:
            if (compIdx == waterCompIdx)
                return 1.0;
            return 0.0;
 
        case oilPhaseIdx:
            if (compIdx == waterCompIdx)
                return 0.0;
            else if (compIdx == oilCompIdx)
                return 1.0 - fluidSystem().convertXoGToxoG(fluidSystem().convertRsToXoG(Rs(), pvtRegionIdx_),
                                                          pvtRegionIdx_);
            else {
                assert(compIdx == gasCompIdx);
                return fluidSystem().convertXoGToxoG(fluidSystem().convertRsToXoG(Rs(), pvtRegionIdx_),
                                                     pvtRegionIdx_);
            }
            break;
 
        case gasPhaseIdx:
            if (compIdx == waterCompIdx)
                return 0.0;
            else if (compIdx == oilCompIdx)
                return fluidSystem().convertXgOToxgO(fluidSystem().convertRvToXgO(Rv(), pvtRegionIdx_),
                                                    pvtRegionIdx_);
            else {
                assert(compIdx == gasCompIdx);
                return 1.0 - fluidSystem().convertXgOToxgO(fluidSystem().convertRvToXgO(Rv(), pvtRegionIdx_),
                                                          pvtRegionIdx_);
            }
            break;
        }
 
        OPM_THROW(std::logic_error, "Invalid phase or component index!");
    }

    OPM_HOST_DEVICE Scalar molarity(unsigned phaseIdx, unsigned compIdx) const
    { return moleFraction(phaseIdx, compIdx)*molarDensity(phaseIdx); }

    OPM_HOST_DEVICE Scalar averageMolarMass(unsigned phaseIdx) const
    {
        Scalar result(0.0);
        for (unsigned compIdx = 0; compIdx < numComponents; ++ compIdx)
            result += fluidSystem().molarMass(compIdx, pvtRegionIdx_)*moleFraction(phaseIdx, compIdx);
        return result;
    }

    OPM_HOST_DEVICE Scalar fugacityCoefficient(unsigned phaseIdx, unsigned compIdx) const
    { return fluidSystem().fugacityCoefficient(*this, phaseIdx, compIdx, pvtRegionIdx_); }

    OPM_HOST_DEVICE Scalar fugacity(unsigned phaseIdx, unsigned compIdx) const
    {
        return
            fugacityCoefficient(phaseIdx, compIdx)
            *moleFraction(phaseIdx, compIdx)
            *pressure(phaseIdx);
    }

    bool phaseIsActive(int phaseIdx) const
    {
        return fluidSystem().phaseIsActive(phaseIdx);
    }

    OPM_HOST_DEVICE const FluidSystem& fluidSystem() const
    {
        if constexpr (fluidSystemIsStatic) {
            static FluidSystem instance;
            return instance;
        } else {
            return **fluidSystemPtr_;
        }
    }
 
private:
    OPM_HOST_DEVICE static unsigned storageToCanonicalPhaseIndex_(unsigned storagePhaseIdx, const FluidSystem& fluidSystem)
    {
        if constexpr (numStoragePhases == 3)
            return storagePhaseIdx;
        else
            return fluidSystem.activeToCanonicalPhaseIdx(storagePhaseIdx);
    }
 
    OPM_HOST_DEVICE static unsigned canonicalToStoragePhaseIndex_(unsigned canonicalPhaseIdx, const FluidSystem& fluidSystem)
    {
        if constexpr (numStoragePhases == 3)
            return canonicalPhaseIdx;
        else
            return fluidSystem.canonicalToActivePhaseIdx(canonicalPhaseIdx);
    }
 
    ConditionalStorage<enableTemperature || enableEnergy, Scalar> temperature_{};
    ConditionalStorage<enableEnergy, std::array<Scalar, numStoragePhases> > enthalpy_{};
    Scalar totalSaturation_{};
    std::array<Scalar, numStoragePhases> pressure_{};
    std::array<Scalar, numStoragePhases> saturation_{};
    std::array<Scalar, numStoragePhases> invB_{};
    std::array<Scalar, numStoragePhases> density_{};
    ConditionalStorage<enableDissolution,Scalar> Rs_{};
    ConditionalStorage<enableDissolution, Scalar> Rv_{};
    ConditionalStorage<enableVapwat,Scalar> Rvw_{};
    ConditionalStorage<enableDissolutionInWater,Scalar> Rsw_{};
    ConditionalStorage<enableBrine, Scalar> saltConcentration_{};
    ConditionalStorage<enableSaltPrecipitation, Scalar> saltSaturation_{};
 
    unsigned short pvtRegionIdx_{};
 
    // If we have a non-static fluid system, we need to store a pointer
    // to it. Otherwise, we do not need to store anything.
    ConditionalStorage<!fluidSystemIsStatic, FluidSystem const*> fluidSystemPtr_;
};
 
} // namespace Opm
 
#endif