EclEpsTwoPhaseLaw.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:
/*
  Copyright (C) 2015 by Andreas Lauser

  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/>.
*/
#ifndef OPM_ECL_EPS_TWO_PHASE_LAW_HPP
#define OPM_ECL_EPS_TWO_PHASE_LAW_HPP

#include "EclEpsTwoPhaseLawParams.hpp"

#include <opm/material/fluidstates/SaturationOverlayFluidState.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/common/Exceptions.hpp>

namespace Opm {
template <class EffLawT,
          class ParamsT = EclEpsTwoPhaseLawParams<EffLawT> >
class EclEpsTwoPhaseLaw : public EffLawT::Traits
{
    typedef EffLawT EffLaw;

public:
    typedef typename EffLaw::Traits Traits;
    typedef ParamsT Params;
    typedef typename EffLaw::Scalar Scalar;

    enum { wettingPhaseIdx = Traits::wettingPhaseIdx };
    enum { nonWettingPhaseIdx = Traits::nonWettingPhaseIdx };

    static const int numPhases = EffLaw::numPhases;
    static_assert(numPhases == 2,
                  "The endpoint scaling applies to the nested twophase laws, not to "
                  "the threephase one!");

    static const bool implementsTwoPhaseApi = true;

    static_assert(EffLaw::implementsTwoPhaseApi,
                  "The material laws put into EclEpsTwoPhaseLaw must implement the "
                  "two-phase material law API!");

    static const bool implementsTwoPhaseSatApi = true;

    static_assert(EffLaw::implementsTwoPhaseSatApi,
                  "The material laws put into EclEpsTwoPhaseLaw must implement the "
                  "two-phase material law saturation API!");

    static const bool isSaturationDependent = true;

    static const bool isPressureDependent = false;

    static const bool isTemperatureDependent = false;

    static const bool isCompositionDependent = false;

    template <class Container, class FluidState>
    static void capillaryPressures(Container& /*values*/, const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The capillaryPressures(fs) method is not yet implemented");
    }

    template <class Container, class FluidState>
    static void relativePermeabilities(Container& /*values*/, const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The pcnw(fs) method is not yet implemented");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcnw(const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The pcnw(fs) method is not yet implemented");
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& SwScaled)
    {
        const Evaluation& SwUnscaled = scaledToUnscaledSatPc(params, SwScaled);
        const Evaluation& pcUnscaled = EffLaw::twoPhaseSatPcnw(params.effectiveLawParams(), SwUnscaled);
        return unscaledToScaledPcnw_(params, pcUnscaled);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatPcnwInv(const Params &params, const Evaluation& pcnwScaled)
    {
        Evaluation pcnwUnscaled = scaledToUnscaledPcnw_(params, pcnwScaled);
        Evaluation SwUnscaled = EffLaw::twoPhaseSatPcnwInv(params.effectiveLawParams(), pcnwUnscaled);
        return unscaledToScaledSatPc(params, SwUnscaled);
    }

    template <class Container, class FluidState>
    static void saturations(Container& /*values*/, const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The saturations(fs) method is not yet implemented");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sw(const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The Sw(fs) method is not yet implemented");
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatSw(const Params& /*params*/, const Evaluation& /*pc*/)
    {
        OPM_THROW(NotImplemented,
                  "The twoPhaseSatSw(pc) method is not yet implemented");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sn(const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The Sn(pc) method is not yet implemented");
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatSn(const Params& /*params*/, const Evaluation& /*pc*/)
    {
        OPM_THROW(NotImplemented,
                  "The twoPhaseSatSn(pc) method is not yet implemented");
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krw(const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The krw(fs) method is not yet implemented");
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& SwScaled)
    {
        const Evaluation& SwUnscaled = scaledToUnscaledSatKrw(params, SwScaled);
        const Evaluation& krwUnscaled = EffLaw::twoPhaseSatKrw(params.effectiveLawParams(), SwUnscaled);
        return unscaledToScaledKrw_(params, krwUnscaled);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatKrwInv(const Params &params, const Evaluation& krwScaled)
    {
        Evaluation krwUnscaled = scaledToUnscaledKrw_(params, krwScaled);
        Evaluation SwUnscaled = EffLaw::twoPhaseSatKrwInv(params.effectiveLawParams(), krwUnscaled);
        return unscaledToScaledSatKrw(params, SwUnscaled);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krn(const Params& /*params*/, const FluidState& /*fluidState*/)
    {
        OPM_THROW(NotImplemented,
                  "The krn(fs) method is not yet implemented");
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& SwScaled)
    {
        const Evaluation& SwUnscaled = scaledToUnscaledSatKrn(params, SwScaled);
        const Evaluation& krnUnscaled = EffLaw::twoPhaseSatKrn(params.effectiveLawParams(), SwUnscaled);
        return unscaledToScaledKrn_(params, krnUnscaled);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatKrnInv(const Params &params, const Evaluation& krnScaled)
    {
        Evaluation krnUnscaled = scaledToUnscaledKrn_(params, krnScaled);
        Evaluation SwUnscaled = EffLaw::twoPhaseSatKrnInv(params.effectiveLawParams(), krnUnscaled);
        return unscaledToScaledSatKrn(params, SwUnscaled);
    }

    template <class Evaluation>
    static Evaluation scaledToUnscaledSatPc(const Params &params, const Evaluation& SwScaled)
    {
        if (!params.config().enableSatScaling())
            return SwScaled;

        // the saturations of capillary pressure are always scaled using two-point
        // scaling
        return scaledToUnscaledSatTwoPoint_(SwScaled,
                                            params.unscaledPoints().saturationPcPoints(),
                                            params.scaledPoints().saturationPcPoints());
    }

    template <class Evaluation>
    static Evaluation unscaledToScaledSatPc(const Params &params, const Evaluation& SwUnscaled)
    {
        if (!params.config().enableSatScaling())
            return SwUnscaled;

        // the saturations of capillary pressure are always scaled using two-point
        // scaling
        return unscaledToScaledSatTwoPoint_(SwUnscaled,
                                            params.unscaledPoints().saturationPcPoints(),
                                            params.scaledPoints().saturationPcPoints());
    }

    template <class Evaluation>
    static Evaluation scaledToUnscaledSatKrw(const Params &params, const Evaluation& SwScaled)
    {
        if (!params.config().enableSatScaling())
            return SwScaled;

        if (params.config().enableThreePointKrSatScaling()) {
            return scaledToUnscaledSatThreePoint_(SwScaled,
                                                  params.unscaledPoints().saturationKrwPoints(),
                                                  params.scaledPoints().saturationKrwPoints());
        }
        else { // two-point relperm saturation scaling
            return scaledToUnscaledSatTwoPoint_(SwScaled,
                                                params.unscaledPoints().saturationKrwPoints(),
                                                params.scaledPoints().saturationKrwPoints());
        }
    }

    template <class Evaluation>
    static Evaluation unscaledToScaledSatKrw(const Params &params, const Evaluation& SwUnscaled)
    {
        if (!params.config().enableSatScaling())
            return SwUnscaled;

        if (params.config().enableThreePointKrSatScaling()) {
            return unscaledToScaledSatThreePoint_(SwUnscaled,
                                                  params.unscaledPoints().saturationKrwPoints(),
                                                  params.scaledPoints().saturationKrwPoints());
        }
        else { // two-point relperm saturation scaling
            return unscaledToScaledSatTwoPoint_(SwUnscaled,
                                                params.unscaledPoints().saturationKrwPoints(),
                                                params.scaledPoints().saturationKrwPoints());
        }
    }

    template <class Evaluation>
    static Evaluation scaledToUnscaledSatKrn(const Params &params, const Evaluation& SwScaled)
    {
        if (!params.config().enableSatScaling())
            return SwScaled;

        if (params.config().enableThreePointKrSatScaling())
            return scaledToUnscaledSatThreePoint_(SwScaled,
                                                  params.unscaledPoints().saturationKrnPoints(),
                                                  params.scaledPoints().saturationKrnPoints());
        else // two-point relperm saturation scaling
            return scaledToUnscaledSatTwoPoint_(SwScaled,
                                                params.unscaledPoints().saturationKrnPoints(),
                                                params.scaledPoints().saturationKrnPoints());
    }


    template <class Evaluation>
    static Evaluation unscaledToScaledSatKrn(const Params &params, const Evaluation& SwUnscaled)
    {
        if (!params.config().enableSatScaling())
            return SwUnscaled;

        if (params.config().enableThreePointKrSatScaling()) {
            return unscaledToScaledSatThreePoint_(SwUnscaled,
                                                  params.unscaledPoints().saturationKrnPoints(),
                                                  params.scaledPoints().saturationKrnPoints());
        }
        else { // two-point relperm saturation scaling
            return unscaledToScaledSatTwoPoint_(SwUnscaled,
                                                params.unscaledPoints().saturationKrnPoints(),
                                                params.scaledPoints().saturationKrnPoints());
        }
    }

private:
    template <class Evaluation, class PointsContainer>
    static Evaluation scaledToUnscaledSatTwoPoint_(const Evaluation& scaledSat,
                                                   const PointsContainer& unscaledSats,
                                                   const PointsContainer& scaledSats)
    {
        return
            unscaledSats[0]
            +
            (scaledSat - scaledSats[0])*((unscaledSats[1] - unscaledSats[0])/(scaledSats[1] - scaledSats[0]));
    }

    template <class Evaluation, class PointsContainer>
    static Evaluation unscaledToScaledSatTwoPoint_(const Evaluation& unscaledSat,
                                                   const PointsContainer& unscaledSats,
                                                   const PointsContainer& scaledSats)
    {
        return
            scaledSats[0]
            +
            (unscaledSat - unscaledSats[0])*((scaledSats[1] - scaledSats[0])/(unscaledSats[1] - unscaledSats[0]));
    }

    template <class Evaluation, class PointsContainer>
    static Evaluation scaledToUnscaledSatThreePoint_(const Evaluation& scaledSat,
                                                     const PointsContainer& unscaledSats,
                                                     const PointsContainer& scaledSats)
    {
        if (unscaledSats[1] >= unscaledSats[2])
            return scaledToUnscaledSatTwoPoint_(scaledSat, unscaledSats, scaledSats);

        if (scaledSat < scaledSats[1])
            return
                unscaledSats[0]
                +
                (scaledSat - scaledSats[0])*((unscaledSats[1] - unscaledSats[0])/(scaledSats[1] - scaledSats[0]));
        else
            return
                unscaledSats[1]
                +
                (scaledSat - scaledSats[1])*((unscaledSats[2] - unscaledSats[1])/(scaledSats[2] - scaledSats[1]));
    }

    template <class Evaluation, class PointsContainer>
    static Evaluation unscaledToScaledSatThreePoint_(const Evaluation& unscaledSat,
                                                     const PointsContainer& unscaledSats,
                                                     const PointsContainer& scaledSats)
    {
        if (unscaledSats[1] >= unscaledSats[2])
            return unscaledToScaledSatTwoPoint_(unscaledSat, unscaledSats, scaledSats);

        if (unscaledSat < unscaledSats[1])
            return
                scaledSats[0]
                +
                (unscaledSat - unscaledSats[0])*((scaledSats[1] - scaledSats[0])/(unscaledSats[1] - unscaledSats[0]));
        else
            return
                scaledSats[1]
                +
                (unscaledSat - unscaledSats[1])*((scaledSats[2] - scaledSats[1])/(unscaledSats[2] - unscaledSats[1]));
    }

    template <class Evaluation>
    static Evaluation unscaledToScaledPcnw_(const Params &params, const Evaluation& unscaledPcnw)
    {
        if (!params.config().enablePcScaling())
            return unscaledPcnw;

        Scalar alpha = params.scaledPoints().maxPcnw()/params.unscaledPoints().maxPcnw();
        return unscaledPcnw*alpha;
    }

    template <class Evaluation>
    static Evaluation scaledToUnscaledPcnw_(const Params &params, const Evaluation& scaledPcnw)
    {
        if (!params.config().enablePcScaling())
            return scaledPcnw;

        Scalar alpha = params.unscaledPoints().maxPcnw()/params.scaledPoints().maxPcnw();
        return scaledPcnw/alpha;
    }

    template <class Evaluation>
    static Evaluation unscaledToScaledKrw_(const Params &params, const Evaluation& unscaledKrw)
    {
        if (!params.config().enableKrwScaling())
            return unscaledKrw;

        // TODO: three point krw y-scaling
        Scalar alpha = params.scaledPoints().maxKrw()/params.unscaledPoints().maxKrw();
        return unscaledKrw*alpha;
    }

    template <class Evaluation>
    static Evaluation scaledToUnscaledKrw_(const Params &params, const Evaluation& scaledKrw)
    {
        if (!params.config().enableKrwScaling())
            return scaledKrw;

        Scalar alpha = params.unscaledPoints().maxKrw()/params.scaledPoints().maxKrw();
        return scaledKrw*alpha;
    }

    template <class Evaluation>
    static Evaluation unscaledToScaledKrn_(const Params &params, const Evaluation& unscaledKrn)
    {
        if (!params.config().enableKrnScaling())
            return unscaledKrn;

        //TODO: three point krn y-scaling
        Scalar alpha = params.scaledPoints().maxKrn()/params.unscaledPoints().maxKrn();
        return unscaledKrn*alpha;
    }

    template <class Evaluation>
    static Evaluation scaledToUnscaledKrn_(const Params &params, const Evaluation& scaledKrn)
    {
        if (!params.config().enableKrnScaling())
            return scaledKrn;

        Scalar alpha = params.unscaledPoints().maxKrn()/params.scaledPoints().maxKrn();
        return scaledKrn*alpha;
    }
};
} // namespace Opm

#endif