ECLPvtCommon.hpp

Go to the documentation of this file.

/*
  Copyright 2017 Statoil ASA.
 
  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 3 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_ECLPVTCOMMON_HEADER_INCLUDED
#define OPM_ECLPVTCOMMON_HEADER_INCLUDED
 
#include <opm/flowdiagnostics/DerivedQuantities.hpp>
#include <opm/utility/ECLPhaseIndex.hpp>
#include <opm/utility/ECLPiecewiseLinearInterpolant.hpp>
#include <opm/utility/ECLPropTable.hpp>
#include <opm/utility/ECLTableInterpolation1D.hpp>
#include <opm/utility/ECLUnitHandling.hpp>
 
#include <algorithm>
#include <cassert>
#include <cstddef>
#include <array>
#include <functional>
#include <initializer_list>
#include <iterator>
#include <memory>
#include <utility>
#include <type_traits>
 
 
namespace Opm {
    class ECLInitFileData;
} // Opm
 
namespace Opm { namespace ECLPVT {
 
    struct ConvertUnits
    {
        using Converter = std::function<double(const double)>;
 
        Converter indep;
 
        std::vector<Converter> column;
    };
 
    struct CreateUnitConverter {
        struct ToSI {
            static ConvertUnits::Converter
            density(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            pressure(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            compressibility(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            disGas(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            vapOil(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvf(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfDerivPress(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfDerivVapOil(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfVisc(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfViscDerivPress(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfViscDerivVapOil(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfGas(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfGasDerivPress(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfGasDerivVapOil(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfGasVisc(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfGasViscDerivPress(const ::Opm::ECLUnits::UnitSystem& usys);
 
            static ConvertUnits::Converter
            recipFvfGasViscDerivVapOil(const ::Opm::ECLUnits::UnitSystem& usys);
        };
    };
 
    enum class RawCurve {
        FVF,
 
        Viscosity,
 
        SaturatedState,
    };
 
    struct PVTGraph {
        std::vector<double> press;
 
        std::vector<double> mixRat;
 
        std::vector<double> value;
    };
 
    template <std::size_t N>
    class DenseVector {
    public:
        explicit DenseVector(const std::array<double, N>& other)
            : x_(other)
        {}
 
        DenseVector& operator+=(const DenseVector& rhs)
        {
            std::transform(std::begin(this->x_),
                           std::end  (this->x_),
                           std::begin(rhs  .x_),
                           std::begin(this->x_),
                           std::plus<double>());
 
            return *this;
        }
 
        DenseVector& operator-=(const DenseVector& rhs)
        {
            std::transform(std::begin(this->x_),
                           std::end  (this->x_),
                           std::begin(rhs  .x_),
                           std::begin(this->x_),
                           std::minus<double>());
 
            return *this;
        }
 
        DenseVector& operator*=(const double rhs)
        {
            std::transform(std::begin(this->x_),
                           std::end  (this->x_),
                           std::begin(this->x_),
               [rhs](const double xi)
            {
                return rhs * xi;
            });
 
            return *this;
        }
 
        DenseVector& operator/=(const double rhs)
        {
            std::transform(std::begin(this->x_),
                           std::end  (this->x_),
                           std::begin(this->x_),
               [rhs](const double xi)
            {
                return xi / rhs;
            });
 
            return *this;
        }
 
        const std::array<double, N>& array() const
        {
            return this->x_;
        }
 
    private:
        std::array<double, N> x_;
    };
 
    template <std::size_t N>
    DenseVector<N> operator/(DenseVector<N> v, const double a)
    {
        return v *= 1.0 / a;
    }
 
    template <std::size_t N>
    DenseVector<N> operator*(const double a, DenseVector<N> v)
    {
        return v *= a;
    }
 
    template <std::size_t N>
    DenseVector<N> operator*(DenseVector<N> v, const double a)
    {
        return v *= a;
    }
 
    template <std::size_t N>
    DenseVector<N> operator+(DenseVector<N> u, const DenseVector<N>& v)
    {
        return u += v;
    }
 
    template <std::size_t N>
    DenseVector<N> operator-(DenseVector<N> u, const DenseVector<N>& v)
    {
        return u -= v;
    }
 
    template <class Extrapolation, bool IsAscendingRange>
    FlowDiagnostics::Graph
    extractRawPVTCurve(const Interp1D::PiecewisePolynomial::Linear<Extrapolation, IsAscendingRange>& interpolant,
                       const RawCurve curve)
    {
        assert ((curve == RawCurve::FVF) ||
                (curve == RawCurve::Viscosity));
 
        const auto colID = (curve == RawCurve::FVF)
            ? std::size_t{0} : std::size_t{1};
 
        auto x = interpolant.independentVariable();
        auto y = interpolant.resultVariable(colID);
 
        assert ((x.size() == y.size()) && "Setup Error");
 
        // Post-process ordinates according to which curve is requested.
        if (curve == RawCurve::FVF) {
            // y == 1/B.  Convert to proper FVF.
            for (auto& yi : y) {
                yi = 1.0 / yi;
            }
        }
        else {
            // y == 1/(B*mu).  Extract viscosity term through the usual
            // conversion formula:
            //
            //    (1 / B) / (1 / (B*mu)).
            const auto b = interpolant.resultVariable(0); // 1/B
 
            assert ((b.size() == y.size()) && "Setup Error");
 
            for (auto n = y.size(), i = 0*n; i < n; ++i) {
                y[i] = b[i] / y[i];
            }
        }
 
        // Graph == pair<vector<double>, vector<double>>
        return FlowDiagnostics::Graph { std::move(x), std::move(y) };
    }
 
    class PVDx
    {
    public:
        using ElemIt = ECLPropTableRawData::ElementIterator;
 
        PVDx(ElemIt               xBegin,
             ElemIt               xEnd,
             const ConvertUnits&  convert,
             std::vector<ElemIt>& colIt);
 
        std::vector<double>
        formationVolumeFactor(const std::vector<double>& p) const;
 
        std::vector<double>
        viscosity(const std::vector<double>& p) const;
 
        PVTGraph getPvtCurve(const RawCurve curve) const;
 
    private:
        using Extrap = ::Opm::Interp1D::PiecewisePolynomial::
            ExtrapolationPolicy::Linearly;
 
        using Backend = ::Opm::Interp1D::PiecewisePolynomial::Linear<Extrap>;
 
        using EvalPt = std::decay<
            decltype(std::declval<Backend>().classifyPoint(0.0))
        >::type;
 
        Backend interp_;
 
        EvalPt getInterpPoint(const double p) const
        {
            return this->interp_.classifyPoint(p);
        }
 
        double fvf_recip(const EvalPt& pt) const
        {
            const auto col = std::size_t{0};
 
            return this->interp_.evaluate(col, pt);
        }
 
        double fvf_mu_recip(const EvalPt& pt) const
        {
            const auto col = std::size_t{1};
 
            return this->interp_.evaluate(col, pt);
        }
 
        template <class EvalDynamicQuant>
        std::vector<double>
        computeQuantity(const std::vector<double>& p,
                        EvalDynamicQuant&&         eval) const
        {
            auto result = std::vector<double>{};
            result.reserve(p.size());
 
            for (const auto& pi : p) {
                result.push_back(eval(this->getInterpPoint(pi)));
            }
 
            return result;
        }
    };
 
    template <class SubtableInterpolant>
    class PVTx
    {
    public:
        PVTx(std::vector<double>              key,
             std::vector<SubtableInterpolant> propInterp)
            : key_       (std::move(key))
            , propInterp_(std::move(propInterp))
        {
            if (this->key_.size() != this->propInterp_.size()) {
                throw std::invalid_argument {
                    "Size of Key Table Does Not Match "
                    "Number of Sub-Table Interpolants"
                };
            }
 
            if (this->key_.size() < 2) {
                throw std::invalid_argument {
                    "Mixing-Dependent Property Evaluator "
                    "Must Have At Least Two Inner Tables"
                };
            }
        }
 
        struct PrimaryKey
        {
            const std::vector<double>& data;
        };
 
        struct InnerVariate
        {
            const std::vector<double>& data;
        };
 
        std::vector<double>
        formationVolumeFactor(const PrimaryKey&   key,
                              const InnerVariate& x) const
        {
            const auto col = std::size_t{0};
 
            return this->computeQuantity(key, x,
                [this, col](const std::size_t     curve,
                            const InnerEvalPoint& pt) -> double
            {   // IFunc: Interpolate 1 / B.
                return this->propInterp_[curve].evaluate(col, pt);
            },
                [](const double recipFvF) -> double
            {
                // OFunc: Convert reciprocal FvF to ordinary FvF.
                return 1.0 / recipFvF;
            });
        }
 
        std::vector<double>
        viscosity(const PrimaryKey& key,
                  const InnerVariate& x) const
        {
            return this->computeQuantity(key, x,
                [this](const std::size_t     curve,
                       const InnerEvalPoint& pt) -> DenseVector<2>
            {
                // IFunc: Interpolate 1/B and 1/(B*mu)
                const auto& I = this->propInterp_[curve];
 
                const auto fvf_recip    = I.evaluate(0, pt);
                const auto fvf_mu_recip = I.evaluate(1, pt);
 
                // { (1 / B) , (1 / (B*mu)) }
                return DenseVector<2>{
                    std::array<double,2>{
                        { fvf_recip, fvf_mu_recip }
                    }
                };
            },
                [](const DenseVector<2>& recipFvFVisc) -> double
            {
                // OFunc: Compute viscosity as
                //   (1 / B) / (1 / B*mu)
 
                const auto& v = recipFvFVisc.array();
 
                return v[0] / v[1];
            });
        }
 
        std::vector<PVTGraph>
        getPvtCurve(const RawCurve curve) const
        {
            if (curve == RawCurve::SaturatedState) {
                return this->saturatedStateCurve();
            }
 
            return this->mainPvtCurve(curve);
        }
 
    private:
        using InnerEvalPoint = typename std::decay<
            decltype(std::declval<SubtableInterpolant>().classifyPoint(0.0))
        >::type;
 
        using OuterInterpPoint =
            ::Opm::Interp1D::PiecewisePolynomial::LocalInterpPoint;
 
        struct InnerInterpPoint {
            InnerEvalPoint left;
            InnerEvalPoint right;
        };
 
        std::vector<double> key_;
        std::vector<SubtableInterpolant> propInterp_;
 
        InnerInterpPoint getInterpPoint(const std::size_t i,
                                        const double      x) const
        {
            assert ((i + 1) < this->propInterp_.size());
 
            return {
                this->propInterp_[i + 0].classifyPoint(x),
                this->propInterp_[i + 1].classifyPoint(x)
            };
        }
 
        template <class Function>
        auto interpolate(Function&&              func,
                         const OuterInterpPoint& outer,
                         const double            x) const
            -> decltype(func(outer.interval, std::declval<InnerEvalPoint>()))
        {
            assert (outer.cat == ::Opm::Interp1D::PointCategory::InRange);
 
            const auto pt =
                this->getInterpPoint(outer.interval, x);
 
            const auto yLeft  = func(outer.interval + 0, pt.left);
            const auto yRight = func(outer.interval + 1, pt.right);
 
            const auto t = outer.t / (this->key_[outer.interval + 1] -
                                      this->key_[outer.interval + 0]);
 
            // t == 0 => yLeft, t == 1 => yRight
            return t*yRight + (1.0 - t)*yLeft;
        }
 
        template <class Function>
        auto extrapLeft(Function&&              func,
                        const OuterInterpPoint& outer,
                        const double            x) const
            -> decltype(func(outer.interval, std::declval<InnerEvalPoint>()))
        {
            assert (outer.cat == ::Opm::Interp1D::PointCategory::LeftOfRange);
            assert (outer.interval == 0*this->key_.size());
 
            const auto pt =
                this->getInterpPoint(outer.interval, x);
 
            const auto yLeft  = func(0, pt.left);
            const auto yRight = func(1, pt.right);
 
            const auto dydk =
                (yRight - yLeft) / (this->key_[1] - this->key_[0]);
 
            return yLeft + outer.t*dydk;
        }
 
        template <class Function>
        auto extrapRight(Function&&              func,
                         const OuterInterpPoint& outer,
                         const double            x) const
            -> decltype(func(outer.interval, std::declval<InnerEvalPoint>()))
        {
            const auto nIntervals = this->key_.size() - 1;
 
            assert (outer.cat == ::Opm::Interp1D::PointCategory::RightOfRange);
            assert (outer.interval == nIntervals);
 
            const auto pt = this->getInterpPoint(nIntervals - 1, x);
 
            const auto yLeft  = func(nIntervals - 1, pt.left);
            const auto yRight = func(nIntervals - 0, pt.right);
 
            const auto dydk =
                (yRight - yLeft) / (this->key_[nIntervals - 0] -
                                    this->key_[nIntervals - 1]);
 
            return yRight + outer.t*dydk;
        }
 
        template <class Function>
        auto evaluate(const double key,
                      const double x,
                      Function&&   func) const
            -> decltype(func(std::declval<OuterInterpPoint>().interval,
                             std::declval<InnerEvalPoint>()))
        {
            const auto outer = ::Opm::Interp1D::PiecewisePolynomial::
                LocalInterpPoint::identify(this->key_, key);
 
            switch (outer.cat) {
            case ::Opm::Interp1D::PointCategory::InRange:
                return this->interpolate(std::forward<Function>(func),
                                         outer, x);
 
            case ::Opm::Interp1D::PointCategory::LeftOfRange:
                return this->extrapLeft(std::forward<Function>(func),
                                        outer, x);
 
            case ::Opm::Interp1D::PointCategory::RightOfRange:
                return this->extrapRight(std::forward<Function>(func),
                                         outer, x);
            }
 
            throw std::logic_error {
                "Outer/Primary Key Cannot Be Classified"
            };
        }
 
        template <class InnerFunction, class OuterFunction>
        std::vector<double>
        computeQuantity(const PrimaryKey&   key,
                        const InnerVariate& x,
                        InnerFunction&&     ifunc,
                        OuterFunction       ofunc) const
        {
            auto result = std::vector<double>{};
 
            const auto nVals = key.data.size();
 
            if (x.data.size() != nVals) {
                throw std::invalid_argument {
                    "Number of Inner Sampling Points Does Not Match "
                    "Number of Outer Sampling Points"
                };
            }
 
            result.reserve(nVals);
 
            for (auto i = 0*nVals; i < nVals; ++i) {
                const auto q =
                    this->evaluate(key.data[i], x.data[i],
                                   std::forward<InnerFunction>(ifunc));
 
                result.push_back(ofunc(q));
            }
 
            return result;
        }
 
        std::vector<PVTGraph>
        mainPvtCurve(const RawCurve curve) const
        {
            // Uses setup appropriate for PVTO.  The PVTG case must be
            // handled in the caller.
 
            auto ret = std::vector<PVTGraph>{};
            ret.reserve(this->propInterp_.size());
 
            auto i = static_cast<std::vector<double>::size_type>(0);
 
            for (const auto& interp : this->propInterp_) {
                auto basic = extractRawPVTCurve(interp, curve);
 
                ret.emplace_back();
                auto& pvtgraph = ret.back();
 
                pvtgraph.mixRat.assign(basic.first.size(), this->key_[i]);
 
                pvtgraph.press = std::move(basic.first);
                pvtgraph.value = std::move(basic.second);
 
                i += 1;
            }
 
            return ret;
        }
 
        std::vector<PVTGraph> saturatedStateCurve() const
        {
            // Uses setup appropriate for PVTO.  The PVTG case must be
            // handled in the caller.
 
            auto curve = PVTGraph{};
 
            curve.press.reserve(this->propInterp_.size());
            curve.mixRat = this->key_;
            curve.value  = this->key_;
 
            for (const auto& interp : this->propInterp_) {
                const auto& yi = interp.independentVariable();
 
                curve.press.push_back(yi[0]);
            }
 
            return { curve };
        }
    };
 
    std::vector<double>
    surfaceMassDensity(const ECLInitFileData& init,
                       const ECLPhaseIndex    phase);
}} // Opm::ECLPVT
 
#endif // OPM_ECLPVTCOMMON_HEADER_INCLUDED