RegularizedVanGenuchten.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) 2008-2013 by Andreas Lauser
  Copyright (C) 2010 by Philipp Nuske
  Copyright (C) 2010 by Bernd Flemisch

  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_REGULARIZED_VAN_GENUCHTEN_HPP
#define OPM_REGULARIZED_VAN_GENUCHTEN_HPP

#include "VanGenuchten.hpp"
#include "RegularizedVanGenuchtenParams.hpp"

#include <opm/material/common/Spline.hpp>

#include <algorithm>

namespace Opm {

template <class TraitsT, class ParamsT = RegularizedVanGenuchtenParams<TraitsT> >
class RegularizedVanGenuchten : public TraitsT
{
    typedef Opm::VanGenuchten<TraitsT, ParamsT> VanGenuchten;

public:
    typedef TraitsT Traits;
    typedef ParamsT Params;

    typedef typename Traits::Scalar Scalar;

    static const int numPhases = Traits::numPhases;
    static_assert(numPhases == 2,
                  "The regularized van Genuchten capillary pressure law only "
                  "applies to the case of two fluid phases");

    static const bool implementsTwoPhaseApi = true;

    static const bool implementsTwoPhaseSatApi = true;

    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 &fs)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;

        values[Traits::wettingPhaseIdx] = 0.0; // reference phase
        values[Traits::nonWettingPhaseIdx] = pcnw<FluidState, Evaluation>(params, fs);
    }

    template <class Container, class FluidState>
    static void saturations(Container &values, const Params &params, const FluidState &fs)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;

        values[Traits::wettingPhaseIdx] = Sw<FluidState, Evaluation>(params, fs);
        values[Traits::nonWettingPhaseIdx] = 1 - values[Traits::wettingPhaseIdx];
    }

    template <class Container, class FluidState>
    static void relativePermeabilities(Container &values, const Params &params, const FluidState &fs)
    {
        typedef typename std::remove_reference<decltype(values[0])>::type Evaluation;

        values[Traits::wettingPhaseIdx] = krw<FluidState, Evaluation>(params, fs);
        values[Traits::nonWettingPhaseIdx] = krn<FluidState, Evaluation>(params, fs);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation pcnw(const Params &params, const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
        return twoPhaseSatPcnw(params, Sw);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatPcnw(const Params &params, const Evaluation& Sw)
    {
        // retrieve the low and the high threshold saturations for the
        // unregularized capillary pressure curve from the parameters
        const Scalar SwThLow = params.pcnwLowSw();
        const Scalar SwThHigh = params.pcnwHighSw();

        // make sure that the capillary pressure observes a derivative
        // != 0 for 'illegal' saturations. This is favourable for the
        // newton solver (if the derivative is calculated numerically)
        // in order to get the saturation moving to the right
        // direction if it temporarily is in an 'illegal' range.
        if (Sw < SwThLow) {
            return params.pcnwLow() + params.pcnwSlopeLow()*(Sw - SwThLow);
        }
        else if (Sw > SwThHigh)
        {
            Scalar yTh = params.pcnwHigh();
            Scalar m1 = (0.0 - yTh)/(1.0 - SwThHigh)*2;

            if (Sw < 1.0) {
                // use spline between threshold Sw and 1.0
                const Spline<Scalar> &sp = params.pcnwHighSpline();

                return sp.eval(Sw);
            }
            else {
                // straight line for Sw > 1.0
                return m1*(Sw - 1.0) + 0.0;
            }
        }

        // if the effective saturation is in an 'reasonable'
        // range, we use the real van genuchten law...
        return VanGenuchten::twoPhaseSatPcnw(params, Sw);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sw(const Params &params, const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& pC =
            FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::nonWettingPhaseIdx))
            - FsToolbox::template toLhs<Evaluation>(fs.pressure(Traits::wettingPhaseIdx));
        return twoPhaseSatSw(params, pC);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatSw(const Params &params, const Evaluation& pC)
    {
        // retrieve the low and the high threshold saturations for the
        // unregularized capillary pressure curve from the parameters
        const Scalar SwThLow = params.pcnwLowSw();
        const Scalar SwThHigh = params.pcnwHighSw();

        // calculate the saturation which corrosponds to the
        // saturation in the non-regularized verision of van
        // Genuchten's law
        if (pC <= 0) {
            // invert straight line for Sw > 1.0
            Scalar m1 = params.pcnwSlopeHigh();
            return pC/m1 + 1.0;
        }

        Evaluation Sw = VanGenuchten::twoPhaseSatSw(params, pC);

        // invert the regularization if necessary
        if (Sw <= SwThLow) {
            // invert the low saturation regularization of pC()
            Scalar pC_SwLow = VanGenuchten::twoPhaseSatPcnw(params, SwThLow);
            return (pC - pC_SwLow)/params.pcnwSlopeLow() + SwThLow;
        }
        else if (SwThHigh < Sw /* && Sw < 1.0*/)
        {
            // invert spline between threshold saturation and 1.0
            const Spline<Scalar>& spline = params.pcnwHighSpline();

            return spline.template intersectInterval<Evaluation>(/*x0=*/SwThHigh, /*x1=*/1.0,
                                                                 /*a=*/0, /*b=*/0, /*c=*/0, /*d=*/pC);
        }

        return Sw;
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation Sn(const Params &params, const FluidState &fs)
    { return 1 - Sw<FluidState, Evaluation>(params, fs); }

    template <class Evaluation>
    static Evaluation twoPhaseSatSn(const Params &params, const Evaluation& pc)
    { return 1 - twoPhaseSatSw(params, pc); }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krw(const Params &params, const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const auto& Sw = FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));
        return twoPhaseSatKrw(params, Sw);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatKrw(const Params &params, const Evaluation& Sw)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        // regularize
        if (Sw <= 0)
            return Toolbox::createConstant(0);
        else if (Sw >= 1)
            return Toolbox::createConstant(1);

        return VanGenuchten::twoPhaseSatKrw(params, Sw);
    }

    template <class FluidState, class Evaluation = typename FluidState::Scalar>
    static Evaluation krn(const Params &params, const FluidState &fs)
    {
        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;

        const Evaluation& Sw =
            1.0 - FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::nonWettingPhaseIdx));
        return twoPhaseSatKrn(params, Sw);
    }

    template <class Evaluation>
    static Evaluation twoPhaseSatKrn(const Params &params, const Evaluation& Sw)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        // regularize
        if (Sw <= 0)
            return Toolbox::createConstant(1);
        else if (Sw >= 1)
            return Toolbox::createConstant(0);

        return VanGenuchten::twoPhaseSatKrn(params, Sw);
    }
};

} // namespace Opm

#endif