EclTwoPhaseMaterial.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_ECL_TWO_PHASE_MATERIAL_HPP
#define OPM_ECL_TWO_PHASE_MATERIAL_HPP

#include <opm/material/fluidmatrixinteractions/EclTwoPhaseMaterialParams.hpp>

#include <opm/common/TimingMacros.hpp>
#include <opm/common/utility/gpuDecorators.hpp>

#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/MathToolbox.hpp>

namespace Opm {

template <class TraitsT,
          class GasOilMaterialLawT,
          class OilWaterMaterialLawT,
          class GasWaterMaterialLawT,
          class ParamsT = EclTwoPhaseMaterialParams<TraitsT,
                                                    typename GasOilMaterialLawT::Params,
                                                    typename OilWaterMaterialLawT::Params,
                                                    typename GasWaterMaterialLawT::Params> >
class EclTwoPhaseMaterial : public TraitsT
{
public:
    using GasOilMaterialLaw = GasOilMaterialLawT;
    using OilWaterMaterialLaw = OilWaterMaterialLawT;
    using GasWaterMaterialLaw = GasWaterMaterialLawT;

    // some safety checks
    static_assert(TraitsT::numPhases == 3,
                  "The number of phases considered by this capillary pressure "
                  "law is always three!");
    static_assert(GasOilMaterialLaw::numPhases == 2,
                  "The number of phases considered by the gas-oil capillary "
                  "pressure law must be two!");
    static_assert(OilWaterMaterialLaw::numPhases == 2,
                  "The number of phases considered by the oil-water capillary "
                  "pressure law must be two!");
    static_assert(GasWaterMaterialLaw::numPhases == 2,
                  "The number of phases considered by the gas-water capillary "
                  "pressure law must be two!");
    static_assert(std::is_same<typename GasOilMaterialLaw::Scalar,
                               typename OilWaterMaterialLaw::Scalar>::value,
                  "The two two-phase capillary pressure laws must use the same "
                  "type of floating point values.");

    using Traits = TraitsT;
    using Params = ParamsT;
    using Scalar = typename Traits::Scalar;

    static constexpr int numPhases = 3;
    static constexpr int waterPhaseIdx = Traits::wettingPhaseIdx;
    static constexpr int oilPhaseIdx = Traits::nonWettingPhaseIdx;
    static constexpr int gasPhaseIdx = Traits::gasPhaseIdx;

    static constexpr bool implementsTwoPhaseApi = false;

    static constexpr bool implementsTwoPhaseSatApi = false;

    static constexpr bool isSaturationDependent = true;

    static constexpr bool isPressureDependent = false;

    static constexpr bool isTemperatureDependent = false;

    static constexpr bool isCompositionDependent = false;

    template <class ContainerT, class FluidState>
    static Scalar relpermOilInOilGasSystem(const Params& /*params*/,
                                           const FluidState& /*fluidState*/) {
        throw std::logic_error {
            "relpermOilInOilGasSystem() is specific to three phases"
                };
    }
    template <class ContainerT, class FluidState>
    static Scalar relpermOilInOilWaterSystem(const Params& /*params*/,
                                                 const FluidState& /*fluidState*/) {
        throw std::logic_error {
                "relpermOilInOilWaterSystem() is specific to three phases"
                    };
    }

    template <class ContainerT, class FluidState, class ...Args>
    OPM_HOST_DEVICE static void capillaryPressures(ContainerT& values,
                                                   const Params& params,
                                                   const FluidState& fluidState)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        using Evaluation = typename std::remove_reference<decltype(values[0])>::type;

        switch (params.approach()) {
        case EclTwoPhaseApproach::GasOil: {
            const Evaluation& So =
                decay<Evaluation>(fluidState.saturation(oilPhaseIdx));

            values[oilPhaseIdx] = 0.0;
            values[gasPhaseIdx] = GasOilMaterialLaw::template twoPhaseSatPcnw<Evaluation, Args...>(params.gasOilParams(), So);
            break;
        }

        case EclTwoPhaseApproach::OilWater: {
            const Evaluation& Sw =
                decay<Evaluation>(fluidState.saturation(waterPhaseIdx));

            values[waterPhaseIdx] = 0.0;
            values[oilPhaseIdx] = OilWaterMaterialLaw::template twoPhaseSatPcnw<Evaluation, Args...>(params.oilWaterParams(), Sw);
            break;
        }

        case EclTwoPhaseApproach::GasWater: {
            const Evaluation& Sw =
                decay<Evaluation>(fluidState.saturation(waterPhaseIdx));

            values[waterPhaseIdx] = 0.0;
            values[gasPhaseIdx] = GasWaterMaterialLaw::template twoPhaseSatPcnw<Evaluation, Args...>(params.gasWaterParams(), Sw);
            break;
        }

        }
    }

    /*
     * Hysteresis parameters for oil-water
     * @see EclHysteresisTwoPhaseLawParams::soMax(...)
     * @see EclHysteresisTwoPhaseLawParams::swMax(...)
     * @see EclHysteresisTwoPhaseLawParams::swMin(...)
     * \param params Parameters
     */
    static void oilWaterHysteresisParams(Scalar& soMax,
                                         Scalar& swMax,
                                         Scalar& swMin,
                                         const Params& params)
    {
        if constexpr (Traits::enableHysteresis) {
            soMax = 1.0 - params.oilWaterParams().krnSwMdc();
            swMax = params.oilWaterParams().krwSwMdc();
            swMin = params.oilWaterParams().pcSwMdc();
            Valgrind::CheckDefined(soMax);
            Valgrind::CheckDefined(swMax);
            Valgrind::CheckDefined(swMin);
        }
    }

    /*
     * Hysteresis parameters for oil-water
     * @see EclHysteresisTwoPhaseLawParams::soMax(...)
     * @see EclHysteresisTwoPhaseLawParams::swMax(...)
     * @see EclHysteresisTwoPhaseLawParams::swMin(...)
     * \param params Parameters
     */
    static void setOilWaterHysteresisParams(const Scalar& soMax,
                                            const Scalar& swMax,
                                            const Scalar& swMin,
                                            Params& params)
    {
        if constexpr (Traits::enableHysteresis) {
            params.oilWaterParams().update(swMin, swMax, 1.0 - soMax);
        }
    }


    /*
     * Hysteresis parameters for gas-oil
     * @see EclHysteresisTwoPhaseLawParams::sgmax(...)
     * @see EclHysteresisTwoPhaseLawParams::shmax(...)
     * @see EclHysteresisTwoPhaseLawParams::somin(...)
     * \param params Parameters
     */
    static void gasOilHysteresisParams(Scalar& sgmax,
                                       Scalar& shmax,
                                       Scalar& somin,
                                       const Params& params)
    {
        if constexpr (Traits::enableHysteresis) {
            sgmax = 1.0 - params.gasOilParams().krnSwMdc();
            shmax = params.gasOilParams().krwSwMdc();
            somin = params.gasOilParams().pcSwMdc();
            Valgrind::CheckDefined(sgmax);
            Valgrind::CheckDefined(shmax);
            Valgrind::CheckDefined(somin);
        }
    }

    /*
     * Hysteresis parameters for gas-oil
     * @see EclHysteresisTwoPhaseLawParams::sgmax(...)
     * @see EclHysteresisTwoPhaseLawParams::shmax(...)
     * \param params Parameters
     */
    static void setGasOilHysteresisParams(const Scalar& sgmax,
                                          const Scalar& shmax,
                                          const Scalar& somin,
                                          Params& params)
    {
        if constexpr (Traits::enableHysteresis) {
            params.gasOilParams().update(somin , shmax, 1.0 - sgmax);
        }
    }

    static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping){
        if(params.approach() == EclTwoPhaseApproach::GasOil)
            return params.gasOilParams().SnTrapped(maximumTrapping);
        if(params.approach() == EclTwoPhaseApproach::GasWater)
            return params.gasWaterParams().SnTrapped(maximumTrapping);
        return 0.0; // oil-water case
    }

    static Scalar strandedGasSaturation(const Params& params, Scalar Sg, Scalar Kg){
        if(params.approach() == EclTwoPhaseApproach::GasOil)
            return params.gasOilParams().SnStranded(Sg, Kg);
        if(params.approach() == EclTwoPhaseApproach::GasWater)
            return params.gasWaterParams().SnStranded(Sg, Kg);
        return 0.0; // oil-water case
    }

    static Scalar trappedOilSaturation(const Params& params, bool maximumTrapping){
        if(params.approach() == EclTwoPhaseApproach::GasOil)
            return params.gasOilParams().SwTrapped();
        if(params.approach() == EclTwoPhaseApproach::OilWater)
            return params.oilWaterParams().SnTrapped(maximumTrapping);
        return 0.0; // gas-water case
    }

    static Scalar trappedWaterSaturation(const Params& params){
        if(params.approach() == EclTwoPhaseApproach::GasWater)
            return params.gasWaterParams().SwTrapped();
        if(params.approach() == EclTwoPhaseApproach::OilWater)
            return params.oilWaterParams().SwTrapped();
        return 0.0; // gas-oil case
    }


    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation pcgn(const Params& /* params */,
                           const FluidState& /* fs */)
    {
        throw std::logic_error("Not implemented: pcgn()");
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation pcnw(const Params& /* params */,
                           const FluidState& /* fs */)
    {
        throw std::logic_error("Not implemented: pcnw()");
    }

    template <class ContainerT, class FluidState>
    static void saturations(ContainerT& /* values */,
                            const Params& /* params */,
                            const FluidState& /* fs */)
    {
        throw std::logic_error("Not implemented: saturations()");
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation Sg(const Params& /* params */,
                         const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: Sg()");
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation Sn(const Params& /* params */,
                         const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: Sn()");
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation Sw(const Params& /* params */,
                         const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: Sw()");
    }

    template <class ContainerT, class FluidState, class ...Args>
    OPM_HOST_DEVICE static void relativePermeabilities(ContainerT& values,
                                                       const Params& params,
                                                       const FluidState& fluidState)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        using Evaluation = typename std::remove_reference<decltype(values[0])>::type;

        switch (params.approach()) {
        case EclTwoPhaseApproach::GasOil: {
            const Evaluation& So =
                decay<Evaluation>(fluidState.saturation(oilPhaseIdx));

            values[oilPhaseIdx] = GasOilMaterialLaw::template twoPhaseSatKrw<Evaluation, Args...>(params.gasOilParams(), So);
            values[gasPhaseIdx] = GasOilMaterialLaw::template twoPhaseSatKrn<Evaluation, Args...>(params.gasOilParams(), So);
            break;
        }

        case EclTwoPhaseApproach::OilWater: {
            const Evaluation& sw =
                decay<Evaluation>(fluidState.saturation(waterPhaseIdx));

            values[waterPhaseIdx] = OilWaterMaterialLaw::template twoPhaseSatKrw<Evaluation, Args...>(params.oilWaterParams(), sw);
            values[oilPhaseIdx] = OilWaterMaterialLaw::template twoPhaseSatKrn<Evaluation, Args...>(params.oilWaterParams(), sw);
            break;
        }

        case EclTwoPhaseApproach::GasWater: {
            const Evaluation& sw =
                decay<Evaluation>(fluidState.saturation(waterPhaseIdx));

            values[waterPhaseIdx] = GasWaterMaterialLaw::template twoPhaseSatKrw<Evaluation, Args...>(params.gasWaterParams(), sw);
            values[gasPhaseIdx] = GasWaterMaterialLaw::template twoPhaseSatKrn<Evaluation, Args...>(params.gasWaterParams(), sw);

            break;
        }
        }
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation krg(const Params& /* params */,
                          const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: krg()");
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation krw(const Params& /* params */,
                          const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: krw()");
    }

    template <class FluidState, class Evaluation = typename FluidState::ValueType>
    static Evaluation krn(const Params& /* params */,
                          const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: krn()");
    }


    template <class FluidState>
    static bool updateHysteresis(Params& params, const FluidState& fluidState)
    {
        if constexpr (Traits::enableHysteresis) {
            OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
            switch (params.approach()) {
            case EclTwoPhaseApproach::GasOil: {
                Scalar So = scalarValue(fluidState.saturation(oilPhaseIdx));
                return params.gasOilParams().update(/*pcSw=*/So, /*krwSw=*/So, /*krnSw=*/So);
            }
            case EclTwoPhaseApproach::OilWater: {
                Scalar sw = scalarValue(fluidState.saturation(waterPhaseIdx));
                return params.oilWaterParams().update(/*pcSw=*/sw, /*krwSw=*/sw, /*krnSw=*/sw);
            }
            case EclTwoPhaseApproach::GasWater: {
                Scalar sw = scalarValue(fluidState.saturation(waterPhaseIdx));
                return params.gasWaterParams().update(/*pcSw=*/sw, /*krwSw=*/sw, /*krnSw=*/sw);
            }
            }
            return false;
        } else {
            return false;
        }
    }
};

} // namespace Opm

#endif