EclMultiplexerMaterial.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_MULTIPLEXER_MATERIAL_HPP
#define OPM_ECL_MULTIPLEXER_MATERIAL_HPP

#include "EclMultiplexerMaterialParams.hpp"
#include "EclDefaultMaterial.hpp"
#include "EclStone1Material.hpp"
#include "EclStone2Material.hpp"
#include "EclTwoPhaseMaterial.hpp"

#include <opm/common/TimingMacros.hpp>
#include <opm/common/ErrorMacros.hpp>

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

#include <algorithm>
#include <stdexcept>

namespace Opm {
#if OPM_IS_INSIDE_DEVICE_FUNCTION
#define OPM_ECL_MULTIPLEXER_MATERIAL_CALL(codeToCall, onePhaseCode)                                                    \
    {                                                                                                                  \
        constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default;                                   \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = DefaultMaterial;                                                                             \
        codeToCall;                                                                                                    \
    }
#else
#define OPM_ECL_MULTIPLEXER_MATERIAL_CALL(codeToCall, onePhaseCode)                                                    \
    switch (params.approach()) {                                                                                       \
    case EclMultiplexerApproach::Stone1: {                                                                             \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Stone1;                   \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = Stone1Material;                                                                              \
        codeToCall;                                                                                                    \
        break;                                                                                                         \
    }                                                                                                                  \
    case EclMultiplexerApproach::Stone2: {                                                                             \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Stone2;                   \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = Stone2Material;                                                                              \
        codeToCall;                                                                                                    \
        break;                                                                                                         \
    }                                                                                                                  \
    case EclMultiplexerApproach::Default: {                                                                            \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default;                  \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = DefaultMaterial;                                                                             \
        codeToCall;                                                                                                    \
        break;                                                                                                         \
    }                                                                                                                  \
    case EclMultiplexerApproach::TwoPhase: {                                                                           \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::TwoPhase;                 \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = TwoPhaseMaterial;                                                                            \
        codeToCall;                                                                                                    \
        break;                                                                                                         \
    }                                                                                                                  \
    case EclMultiplexerApproach::OnePhase: {                                                                           \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::OnePhase;                 \
        onePhaseCode;                                                                                                  \
        break;                                                                                                         \
    }                                                                                                                  \
    }
#endif

// The static_assert does not compile with gcc 12 and earlier (or nvcc) when placed in the multiplexer below.
#if (defined(__GNUC__) && (__GNUC__ < 13)) || defined(__CUDACC__)
    #define STATIC_ASSERT_ECL_MULTIPLEXER_UNLESS_GCC_LT_13 throw std::logic_error("Unhandled EclMultiplexerApproach")
#else
    #define STATIC_ASSERT_ECL_MULTIPLEXER_UNLESS_GCC_LT_13 static_assert(false, "Unhandled EclMultiplexerApproach")
#endif

#if OPM_IS_INSIDE_DEVICE_FUNCTION
#define OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(codeToCall, onePhaseCode)                                        \
    if constexpr (Head::approach == EclMultiplexerApproach::Default) {                                                 \
        constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default;                                   \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = DefaultMaterial;                                                                             \
        codeToCall;                                                                                                    \
    }

#else
#define OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(codeToCall, onePhaseCode)                                        \
    if constexpr (Head::approach == EclMultiplexerApproach::Stone1) {                                                  \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Stone1;                   \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = Stone1Material;                                                                              \
        codeToCall;                                                                                                    \
    } else if constexpr (Head::approach == EclMultiplexerApproach::Stone2) {                                           \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Stone2;                   \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = Stone2Material;                                                                              \
        codeToCall;                                                                                                    \
    } else if constexpr (Head::approach == EclMultiplexerApproach::Default) {                                          \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::Default;                  \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = DefaultMaterial;                                                                             \
        codeToCall;                                                                                                    \
    } else if constexpr (Head::approach == EclMultiplexerApproach::TwoPhase) {                                         \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::TwoPhase;                 \
        auto& realParams = params.template getRealParams<approach>();                                                  \
        using ActualLaw = TwoPhaseMaterial;                                                                            \
        codeToCall;                                                                                                    \
    } else if constexpr (Head::approach == EclMultiplexerApproach::OnePhase) {                                         \
        [[maybe_unused]] constexpr EclMultiplexerApproach approach = EclMultiplexerApproach::OnePhase;                 \
        onePhaseCode;                                                                                                  \
    } else {                                                                                                           \
        STATIC_ASSERT_ECL_MULTIPLEXER_UNLESS_GCC_LT_13;                                                                \
    }
#endif
// Pass this for the onePhaseCode argument if nothing is to be done.
inline void doNothing() { };

template <class TraitsT,
          class GasOilMaterialLawT,
          class OilWaterMaterialLawT,
          class GasWaterMaterialLawT,
          class ParamsT = EclMultiplexerMaterialParams<TraitsT,
                                                       GasOilMaterialLawT,
                                                       OilWaterMaterialLawT,
                                                       GasWaterMaterialLawT> >
class EclMultiplexerMaterial : public TraitsT
{
public:
    using GasOilMaterialLaw = GasOilMaterialLawT;
    using OilWaterMaterialLaw = OilWaterMaterialLawT;
    using GasWaterMaterialLaw = GasWaterMaterialLawT;

    using Stone1Material = EclStone1Material<TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw>;
    using Stone2Material = EclStone2Material<TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw>;
    using DefaultMaterial = EclDefaultMaterial<TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw>;
    using TwoPhaseMaterial = EclTwoPhaseMaterial<TraitsT, GasOilMaterialLaw, OilWaterMaterialLaw, GasWaterMaterialLaw>;

    // 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, class ...Args>
    OPM_HOST_DEVICE static void capillaryPressures(ContainerT& values,
                                                   const Params& params,
                                                   const FluidState& fluidState)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        if constexpr (FrontIsEclMultiplexerDispatchV<Args...>) {
            capillaryPressuresT<ContainerT, FluidState, Args...>(values, params, fluidState);
            return;
        }
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(ActualLaw::capillaryPressures(values, realParams, fluidState),
                                          values[0] = 0.0);
    }


    template <class ContainerT, class FluidState, class Head, class ...Args>
    static void capillaryPressuresT(ContainerT& values,
                                    const Params& params,
                                    const FluidState& fluidState)
    {
#define OPM_LOCAL_TEMPLATE_ARGS ContainerT, FluidState, Args...
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(
            ActualLaw::template capillaryPressures<OPM_LOCAL_TEMPLATE_ARGS>(values, realParams, fluidState),
            values[0] = 0.0
        );
#undef OPM_LOCAL_TEMPLATE_ARGS
    }

    /*
     * 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)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(ActualLaw::oilWaterHysteresisParams(soMax, swMax, swMin, realParams),
                                          doNothing());
    }

    /*
     * 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)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(ActualLaw::setOilWaterHysteresisParams(soMax, swMax, swMin, realParams),
                                          doNothing());
    }

    /*
     * 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)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(ActualLaw::gasOilHysteresisParams(sgmax, shmax, somin, realParams),
                                          doNothing());
    }

    static Scalar trappedGasSaturation(const Params& params, bool maximumTrapping)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(return ActualLaw::trappedGasSaturation(realParams, maximumTrapping),
                                          return 0.0);
        return 0.0;
    }

    static Scalar strandedGasSaturation(const Params& params, Scalar Sg, Scalar Kg)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(return ActualLaw::strandedGasSaturation(realParams, Sg, Kg),
                                          return 0.0);
        return 0.0;
    }

    static Scalar trappedOilSaturation(const Params& params, bool maximumTrapping)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(return ActualLaw::trappedOilSaturation(realParams, maximumTrapping),
                                          return 0.0);
        return 0.0;
    }

    static Scalar trappedWaterSaturation(const Params& params)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(return ActualLaw::trappedWaterSaturation(realParams),
                                          return 0.0);
        return 0.0;
    }
    /*
     * Hysteresis parameters for gas-oil
     * @see EclHysteresisTwoPhaseLawParams::sgmax(...)
     * @see EclHysteresisTwoPhaseLawParams::shmax(...)
     * @see EclHysteresisTwoPhaseLawParams::somin(...)
     * \param params Parameters
     */
    static void setGasOilHysteresisParams(const Scalar& sgmax,
                                          const Scalar& shmax,
                                          const Scalar& somin,
                                          Params& params)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(ActualLaw::setGasOilHysteresisParams(sgmax, shmax, somin, realParams),
                                          doNothing());
    }

    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);
        if constexpr (FrontIsEclMultiplexerDispatchV<Args...>) {
            relativePermeabilitiesT<ContainerT, FluidState, Args...>(values, params, fluidState);
            return;
        }
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(ActualLaw::relativePermeabilities(values, realParams, fluidState),
                                          values[0] = 1.0);
    }

    template <class ContainerT, class FluidState, class Head, class ...Args>
    static void relativePermeabilitiesT(ContainerT& values,
                                        const Params& params,
                                        const FluidState& fluidState)
    {
#define OPM_LOCAL_TEMPLATE_ARGS ContainerT, FluidState, Args...
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL_COMPILETIME(
            ActualLaw::template relativePermeabilities<OPM_LOCAL_TEMPLATE_ARGS>(values, realParams, fluidState),
            values[0] = 1.0
        );
#undef OPM_LOCAL_TEMPLATE_ARGS
    }




    template <class Evaluation, class FluidState>
    static Evaluation relpermOilInOilGasSystem(const Params& params,
                                               const FluidState& fluidState)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        switch (params.approach()) {
        case EclMultiplexerApproach::Stone1:
            return Stone1Material::template relpermOilInOilGasSystem<Evaluation>
                (params.template getRealParams<EclMultiplexerApproach::Stone1>(),
                 fluidState);

        case EclMultiplexerApproach::Stone2:
            return Stone2Material::template relpermOilInOilGasSystem<Evaluation>
                (params.template getRealParams<EclMultiplexerApproach::Stone2>(),
                 fluidState);

        case EclMultiplexerApproach::Default:
            return DefaultMaterial::template relpermOilInOilGasSystem<Evaluation>
                (params.template getRealParams<EclMultiplexerApproach::Default>(),
                 fluidState);

        default:
            throw std::logic_error {
                "relpermOilInOilGasSystem() is specific to three phases"
            };
        }
    }

    template <class Evaluation, class FluidState>
    static Evaluation relpermOilInOilWaterSystem(const Params& params,
                                                 const FluidState& fluidState)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        switch (params.approach()) {
        case EclMultiplexerApproach::Stone1:
            return Stone1Material::template relpermOilInOilWaterSystem<Evaluation>
                (params.template getRealParams<EclMultiplexerApproach::Stone1>(),
                 fluidState);

        case EclMultiplexerApproach::Stone2:
            return Stone2Material::template relpermOilInOilWaterSystem<Evaluation>
                (params.template getRealParams<EclMultiplexerApproach::Stone2>(),
                 fluidState);

        case EclMultiplexerApproach::Default:
            return DefaultMaterial::template relpermOilInOilWaterSystem<Evaluation>
                (params.template getRealParams<EclMultiplexerApproach::Default>(),
                 fluidState);

        default:
            throw std::logic_error {
                "relpermOilInOilWaterSystem() is specific to three phases"
            };
        }
    }

    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)
    {
        OPM_TIMEFUNCTION_LOCAL(Subsystem::SatProps);
        OPM_ECL_MULTIPLEXER_MATERIAL_CALL(return ActualLaw::updateHysteresis(realParams, fluidState),
                                          return false);
        return false;
    }
};

} // namespace Opm

#endif