BrooksCorey.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) 2009-2013 by Andreas Lauser
  Copyright (C) 2011 by Holger Class
  Copyright (C) 2010 by Philipp Nuske

  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_BROOKS_COREY_HPP
#define OPM_BROOKS_COREY_HPP

#include "BrooksCoreyParams.hpp"

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

#include <algorithm>
#include <cassert>
#include <cmath>

namespace Opm {
template <class TraitsT, class ParamsT = BrooksCoreyParams<TraitsT> >
class BrooksCorey : public TraitsT
{
public:
    typedef TraitsT Traits;
    typedef ParamsT Params;
    typedef typename Traits::Scalar Scalar;

    static const int numPhases = Traits::numPhases;
    static_assert(numPhases == 2,
                  "The Brooks-Corey 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;

    static_assert(Traits::numPhases == 2,
                  "The number of fluid phases must be two if you want to use "
                  "this material law!");

    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 Evaluation& Sw =
            FsToolbox::template toLhs<Evaluation>(fs.saturation(Traits::wettingPhaseIdx));

        assert(0 <= Sw && Sw <= 1);

        return twoPhaseSatPcnw(params, Sw);
    }

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

        assert(0 <= Sw && Sw <= 1);

        return params.entryPressure()*Toolbox::pow(Sw, -1/params.lambda());
    }

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

        assert(pcnw > 0.0);

        return Toolbox::pow(params.entryPressure()/pcnw, -params.lambda());
    }

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

        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)
    {
        typedef MathToolbox<Evaluation> Toolbox;

        assert(pc > 0); // if we don't assume that, std::pow will screw up!

        return Toolbox::pow(pc/params.entryPressure(), -params.lambda());
    }

    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;

        assert(0 <= Sw && Sw <= 1);

        return Toolbox::pow(Sw, 2.0/params.lambda() + 3);
    }

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

        return Toolbox::pow(krw, 1.0/(2.0/params.lambda() + 3));
    }

    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;

        assert(0 <= Sw && Sw <= 1);

        Scalar exponent = 2.0/params.lambda() + 1;
        const Evaluation Sn = 1. - Sw;
        return Sn*Sn*(1. - Toolbox::pow(Sw, exponent));
    }

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

        // since inverting the formula for krn is hard to do analytically, we use the
        // Newton-Raphson method
        Evaluation Sw = Toolbox::createConstant(0.5);
        Scalar eps = 1e-10;
        for (int i = 0; i < 20; ++i) {
            Evaluation f = krn - twoPhaseSatKrn(params, Sw);
            Evaluation fStar = krn - twoPhaseSatKrn(params, Sw + eps);
            Evaluation fPrime = (fStar - f)/eps;

            Evaluation delta = f/fPrime;
            Sw -= delta;
            if (Sw < 0)
                Sw = Toolbox::createConstant(0.0);
            if (Toolbox::abs(delta) < 1e-10)
                return Sw;
        }

        OPM_THROW(NumericalProblem,
                  "Couldn't invert the Brooks-Corey non-wetting phase"
                  " relperm within 20 iterations");
    }
};
} // namespace Opm

#endif // BROOKS_COREY_HPP