EclStone1Material.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_STONE1_MATERIAL_HPP
#define OPM_ECL_STONE1_MATERIAL_HPP
 
#include "EclStone1MaterialParams.hpp"
 
#include <opm/material/common/Valgrind.hpp>
#include <opm/material/common/MathToolbox.hpp>
 
#include <algorithm>
#include <cmath>
#include <stdexcept>
#include <type_traits>
 
namespace Opm {
 
template <class TraitsT,
          class GasOilMaterialLawT,
          class OilWaterMaterialLawT,
          class ParamsT = EclStone1MaterialParams<TraitsT, GasOilMaterialLawT, OilWaterMaterialLawT> >
class EclStone1Material : public TraitsT
{
public:
    using GasOilMaterialLaw = GasOilMaterialLawT;
    using OilWaterMaterialLaw = OilWaterMaterialLawT;
 
    // 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(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.");
 
    static_assert(GasOilMaterialLaw::implementsTwoPhaseSatApi,
                  "The gas-oil material law must implement the two-phase saturation "
                  "only API to for the default Ecl capillary pressure law!");
    static_assert(OilWaterMaterialLaw::implementsTwoPhaseSatApi,
                  "The oil-water material law must implement the two-phase saturation "
                  "only API to for the default Ecl capillary pressure law!");
 
    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 void capillaryPressures(ContainerT& values,
                                   const Params& params,
                                   const FluidState& state)
    {
        using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
        values[gasPhaseIdx] = pcgn<FluidState, Evaluation>(params, state);
        values[oilPhaseIdx] = 0;
        values[waterPhaseIdx] = - pcnw<FluidState, Evaluation>(params, state);
        Valgrind::CheckDefined(values[gasPhaseIdx]);
        Valgrind::CheckDefined(values[oilPhaseIdx]);
        Valgrind::CheckDefined(values[waterPhaseIdx]);
    }
 
    /*
     * Hysteresis parameters for oil-water
     * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
     * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
     * \param params Parameters
     */
    static void oilWaterHysteresisParams(Scalar& pcSwMdc,
                                         Scalar& krnSwMdc,
                                         const Params& params)
    {
        pcSwMdc = params.oilWaterParams().pcSwMdc();
        krnSwMdc = params.oilWaterParams().krnSwMdc();
 
        Valgrind::CheckDefined(pcSwMdc);
        Valgrind::CheckDefined(krnSwMdc);
    }
 
    /*
     * Hysteresis parameters for oil-water
     * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
     * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
     * \param params Parameters
     */
    static void setOilWaterHysteresisParams(const Scalar& pcSwMdc,
                                            const Scalar& krnSwMdc,
                                            Params& params)
    {
        constexpr const double krwSw = 2.0; //Should not be used
        params.oilWaterParams().update(pcSwMdc, krwSw, krnSwMdc);
    }
 
    /*
     * Hysteresis parameters for gas-oil
     * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
     * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
     * \param params Parameters
     */
    static void gasOilHysteresisParams(Scalar& pcSwMdc,
                                       Scalar& krnSwMdc,
                                       const Params& params)
    {
        const auto Swco = params.Swl();
 
        // Pretend oil saturation is 'Swco' larger than it really is in
        // order to infer correct maximum Sg values in output layer.
        pcSwMdc = std::min(params.gasOilParams().pcSwMdc() + Swco, Scalar{2.0});
        krnSwMdc = std::min(params.gasOilParams().krnSwMdc() + Swco, Scalar{2.0});
 
        Valgrind::CheckDefined(pcSwMdc);
        Valgrind::CheckDefined(krnSwMdc);
    }
 
    /*
     * Hysteresis parameters for gas-oil
     * @see EclHysteresisTwoPhaseLawParams::pcSwMdc(...)
     * @see EclHysteresisTwoPhaseLawParams::krnSwMdc(...)
     * \param params Parameters
     */
    static void setGasOilHysteresisParams(const Scalar& pcSwMdc,
                                          const Scalar& krnSwMdc,
                                          Params& params)
    {
        // Maximum attainable oil saturation is 1-SWL
        const auto Swco = params.Swl();
        constexpr const double krwSw = 2.0; //Should not be used
        params.gasOilParams().update(pcSwMdc - Swco, krwSw, krnSwMdc - Swco);
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcgn(const Params& params,
                           const FluidState& fs)
    {
        // Maximum attainable oil saturation is 1-SWL
        const auto Sw = 1.0 - params.Swl() - decay<Evaluation>(fs.saturation(gasPhaseIdx));
        return GasOilMaterialLaw::twoPhaseSatPcnw(params.gasOilParams(), Sw);
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcnw(const Params& params,
                           const FluidState& fs)
    {
        const auto Sw = decay<Evaluation>(fs.saturation(waterPhaseIdx));
        Valgrind::CheckDefined(Sw);
 
        const auto result = OilWaterMaterialLaw::twoPhaseSatPcnw(params.oilWaterParams(), Sw);
        Valgrind::CheckDefined(result);
 
        return result;
    }
 
    template <class ContainerT, class FluidState>
    static void saturations(ContainerT& /* values */,
                            const Params& /* params */,
                            const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: saturations()");
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sg(const Params& /* params */,
                         const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: Sg()");
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sn(const Params& /* params */,
                         const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: Sn()");
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sw(const Params& /* params */,
                         const FluidState& /* fluidState */)
    {
        throw std::logic_error("Not implemented: Sw()");
    }
 
    template <class ContainerT, class FluidState>
    static void relativePermeabilities(ContainerT& values,
                                       const Params& params,
                                       const FluidState& fluidState)
    {
        using Evaluation = typename std::remove_reference<decltype(values[0])>::type;
 
        values[waterPhaseIdx] = krw<FluidState, Evaluation>(params, fluidState);
        values[oilPhaseIdx] = krn<FluidState, Evaluation>(params, fluidState);
        values[gasPhaseIdx] = krg<FluidState, Evaluation>(params, fluidState);
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krg(const Params& params,
                          const FluidState& fluidState)
    {
        // Maximum attainable oil saturation is 1-SWL,
        const Evaluation Sw = 1 - params.Swl() - decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
        return GasOilMaterialLaw::twoPhaseSatKrn(params.gasOilParams(), Sw);
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krw(const Params& params,
                          const FluidState& fluidState)
    {
        const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
        return OilWaterMaterialLaw::twoPhaseSatKrw(params.oilWaterParams(), Sw);
    }
 
    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krn(const Params& params,
                          const FluidState& fluidState)
    {
        // the Eclipse docu is inconsistent with naming the variable of connate water: In
        // some places the connate water saturation is represented by "Swl", in others
        // "Swco" is used.
        const Scalar Swco = params.Swl();
 
        // oil relperm at connate water saturations (with Sg=0)
        const Scalar krocw = params.krocw();
 
        const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
        const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
 
        const Evaluation kro_ow = relpermOilInOilWaterSystem<Evaluation>(params, fluidState);
        const Evaluation kro_go = relpermOilInOilGasSystem<Evaluation>(params, fluidState);
 
        Evaluation beta;
        if (Sw <= Swco)
            beta = 1.0;
        else {
            // there seems to be an error in the ECL documentation: using the approach to
            // the scaled saturations as described there leads to significant deviations
            // from the results produced by Eclipse 100.
            const Evaluation SSw = (Sw - Swco)/(1.0 - Swco);
            const Evaluation SSg = Sg/(1.0 - Swco);
            const Evaluation SSo = 1.0 - SSw - SSg;
 
            if (SSw >= 1.0 || SSg >= 1.0)
                beta = 1.0;
            else
                beta = pow( SSo/((1 - SSw)*(1 - SSg)), params.eta());
        }
 
        return max(0.0, min(1.0, beta*kro_ow*kro_go/krocw));
    }
 
    template <class Evaluation, class FluidState>
    static Evaluation relpermOilInOilGasSystem(const Params& params,
                                               const FluidState& fluidState)
    {
        const Evaluation Sg = decay<Evaluation>(fluidState.saturation(gasPhaseIdx));
 
        return GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - params.Swl());
    }
 
    template <class Evaluation, class FluidState>
    static Evaluation relpermOilInOilWaterSystem(const Params& params,
                                                 const FluidState& fluidState)
    {
        const Evaluation Sw = decay<Evaluation>(fluidState.saturation(waterPhaseIdx));
 
        return OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw);
    }
 
    template <class FluidState>
    static void updateHysteresis(Params& params, const FluidState& fluidState)
    {
        const Scalar Swco = params.Swl();
        const Scalar Sw = scalarValue(fluidState.saturation(waterPhaseIdx));
        const Scalar Sg = scalarValue(fluidState.saturation(gasPhaseIdx));
 
        params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/Sw, /*krnSw=*/Sw);
        params.gasOilParams().update(/*pcSw=*/  1.0 - Swco - Sg,
                                     /*krwSw=*/ 1.0 - Swco - Sg,
                                     /*krnSw=*/ 1.0 - Swco - Sg);
    }
};
 
} // namespace Opm
 
#endif