RateConverter.hpp

Go to the documentation of this file.

/*
  Copyright 2014, 2015 SINTEF ICT, Applied Mathematics.
  Copyright 2014, 2015 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_RATECONVERTER_HPP_HEADER_INCLUDED
#define OPM_RATECONVERTER_HPP_HEADER_INCLUDED

#include <opm/autodiff/BlackoilPropsAdInterface.hpp>

#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/RegionMapping.hpp>

#include <Eigen/Core>

#include <algorithm>
#include <cmath>
#include <memory>
#include <stdexcept>
#include <type_traits>
#include <unordered_map>
#include <utility>
#include <vector>

namespace Opm {
    namespace RateConverter {
        namespace Details {
            namespace Select {
                template <class RegionID, bool>
                struct RegionIDParameter
                {
                    using type =
                        typename std::remove_reference<RegionID>::type &;
                };

                template <class RegionID>
                struct RegionIDParameter<RegionID, true>
                {
                    using type = RegionID;
                };
            } // Select

            template <typename RegionId, class Attributes>
            class RegionAttributes
            {
            public:
                using RegionID =
                    typename Select::RegionIDParameter
                    <RegionId, std::is_integral<RegionId>::value>::type;

                template <class RMap>
                RegionAttributes(const RMap&       rmap,
                                 const Attributes& attr)
                {
                    using VT = typename AttributeMap::value_type;

                    for (const auto& r : rmap.activeRegions()) {
                        auto v = std::unique_ptr<Value>(new Value(attr));

                        const auto stat = attr_.insert(VT(r, std::move(v)));

                        if (stat.second) {
                            // New value inserted.
                            const auto& cells = rmap.cells(r);

                            assert (! cells.empty());

                            // Region's representative cell.
                            stat.first->second->cell_ = cells[0];
                        }
                    }
                }

                int cell(const RegionID reg) const
                {
                    return this->find(reg).cell_;
                }

                const Attributes& attributes(const RegionID reg) const
                {
                    return this->find(reg).attr_;
                }

                Attributes& attributes(const RegionID reg)
                {
                    return this->find(reg).attr_;
                }

            private:
                struct Value {
                    Value(const Attributes& attr)
                        : attr_(attr)
                        , cell_(-1)
                    {}

                    Attributes attr_;
                    int        cell_;
                };

                using ID =
                    typename std::remove_reference<RegionId>::type;

                using AttributeMap =
                    std::unordered_map<ID, std::unique_ptr<Value>>;

                AttributeMap attr_;

                const Value& find(const RegionID reg) const
                {
                    const auto& i = attr_.find(reg);

                    if (i == attr_.end()) {
                        throw std::invalid_argument("Unknown region ID");
                    }

                    return *i->second;
                }

                Value& find(const RegionID reg)
                {
                    const auto& i = attr_.find(reg);

                    if (i == attr_.end()) {
                        throw std::invalid_argument("Unknown region ID");
                    }

                    return *i->second;
                }
            };

            namespace PhaseUsed {
                inline bool
                water(const PhaseUsage& pu)
                {
                    return pu.phase_used[ BlackoilPhases::Aqua ] != 0;
                }

                inline bool
                oil(const PhaseUsage& pu)
                {
                    return pu.phase_used[ BlackoilPhases::Liquid ] != 0;
                }

                inline bool
                gas(const PhaseUsage& pu)
                {
                    return pu.phase_used[ BlackoilPhases::Vapour ] != 0;
                }
            } // namespace PhaseUsed

            namespace PhasePos {
                inline int
                water(const PhaseUsage& pu)
                {
                    int p = -1;

                    if (PhaseUsed::water(pu)) {
                        p = pu.phase_pos[ BlackoilPhases::Aqua ];
                    }

                    return p;
                }

                inline int
                oil(const PhaseUsage& pu)
                {
                    int p = -1;

                    if (PhaseUsed::oil(pu)) {
                        p = pu.phase_pos[ BlackoilPhases::Liquid ];
                    }

                    return p;
                }

                inline int
                gas(const PhaseUsage& pu)
                {
                    int p = -1;

                    if (PhaseUsed::gas(pu)) {
                        p = pu.phase_pos[ BlackoilPhases::Vapour ];
                    }

                    return p;
                }
            } // namespace PhasePos
        } // namespace Details

        template <class Property, class Region>
        class SurfaceToReservoirVoidage {
        public:
            SurfaceToReservoirVoidage(const Property& props,
                                      const Region&   region)
                : props_(props)
                , rmap_ (region)
                , attr_ (rmap_, Attributes(props_.numPhases()))
            {}

            void
            defineState(const BlackoilState& state)
            {
                calcAverages(state);
                calcRmax();
            }

            typedef typename RegionMapping<Region>::RegionId RegionId;

            template <class Input,
                      class Coeff>
            void
            calcCoeff(const Input& in, const RegionId r, Coeff& coeff) const
            {
                using V = typename Property::V;

                const auto& pu = props_.phaseUsage();
                const auto& ra = attr_.attributes(r);

                const auto  p  = this->constant(ra.pressure);
                const auto  T  = this->constant(ra.temperature);
                const auto  c  = this->getRegCell(r);

                const int   iw = Details::PhasePos::water(pu);
                const int   io = Details::PhasePos::oil  (pu);
                const int   ig = Details::PhasePos::gas  (pu);

                std::fill(& coeff[0], & coeff[0] + props_.numPhases(), 0.0);

                if (Details::PhaseUsed::water(pu)) {
                    // q[w]_r = q[w]_s / bw

                    const V bw = props_.bWat(p, T, c).value();

                    coeff[iw] = 1.0 / bw(0);
                }

                const Miscibility& m = calcMiscibility(in, r);

                // Determinant of 'R' matrix
                const double detR = 1.0 - (m.rs(0) * m.rv(0));

                if (Details::PhaseUsed::oil(pu)) {
                    // q[o]_r = 1/(bo * (1 - rs*rv)) * (q[o]_s - rv*q[g]_s)

                    const auto   rs  = this->constant(m.rs);
                    const V      bo  = props_.bOil(p, T, rs, m.cond, c).value();
                    const double den = bo(0) * detR;

                    coeff[io] += 1.0 / den;

                    if (Details::PhaseUsed::gas(pu)) {
                        coeff[ig] -= m.rv(0) / den;
                    }
                }

                if (Details::PhaseUsed::gas(pu)) {
                    // q[g]_r = 1/(bg * (1 - rs*rv)) * (q[g]_s - rs*q[o]_s)

                    const auto   rv  = this->constant(m.rv);
                    const V      bg  = props_.bGas(p, T, rv, m.cond, c).value();
                    const double den = bg(0) * detR;

                    coeff[ig] += 1.0 / den;

                    if (Details::PhaseUsed::oil(pu)) {
                        coeff[io] -= m.rs(0) / den;
                    }
                }
            }

        private:
            const Property&  props_;

            const RegionMapping<Region> rmap_;

            struct Attributes {
                Attributes(const int np)
                    : pressure   (0.0)
                    , temperature(0.0)
                    , Rmax       (Eigen::ArrayXd::Zero(np, 1))
                {}

                double         pressure;
                double         temperature;
                Eigen::ArrayXd Rmax;
            };

            Details::RegionAttributes<RegionId, Attributes> attr_;

            struct Miscibility {
                Miscibility()
                    : rs  (1)
                    , rv  (1)
                    , cond(1)
                {
                    rs << 0.0;
                    rv << 0.0;
                }

                typename Property::V rs;

                typename Property::V rv;

                std::vector<PhasePresence> cond;
            };

            void
            calcAverages(const BlackoilState& state)
            {
                const auto& press = state.pressure();
                const auto& temp  = state.temperature();

                for (const auto& reg : rmap_.activeRegions()) {
                    auto& ra = attr_.attributes(reg);
                    auto& p  = ra.pressure;
                    auto& T  = ra.temperature;

                    std::size_t n = 0;
                    p = T = 0.0;
                    for (const auto& cell : rmap_.cells(reg)) {
                        p += press[ cell ];
                        T += temp [ cell ];

                        n += 1;
                    }

                    p /= n;
                    T /= n;
                }
            }

            void
            calcRmax()
            {
                const PhaseUsage& pu = props_.phaseUsage();

                if (Details::PhaseUsed::oil(pu) &&
                    Details::PhaseUsed::gas(pu))
                {
                    const Eigen::ArrayXd::Index
                        io = Details::PhasePos::oil(pu),
                        ig = Details::PhasePos::gas(pu);

                    // Note: Intentionally does not take capillary
                    // pressure into account.  This facility uses the
                    // average *hydrocarbon* pressure rather than
                    // average phase pressure.

                    for (const auto& reg : rmap_.activeRegions()) {
                        auto& ra = attr_.attributes(reg);

                        const auto c = this->getRegCell(reg);
                        const auto p = this->constant(ra.pressure);
                        const auto T = this->constant(ra.temperature);

                        auto& Rmax = ra.Rmax;
                        Rmax.row(io) = props_.rsSat(p, T, c).value();
                        Rmax.row(ig) = props_.rvSat(p, T, c).value();
                    }
                }
            }

            template <class Input>
            Miscibility
            calcMiscibility(const Input& in, const RegionId r) const
            {
                const auto& pu   = props_.phaseUsage();
                const auto& attr = attr_.attributes(r);

                const int io = Details::PhasePos::oil(pu);
                const int ig = Details::PhasePos::gas(pu);

                Miscibility m;
                PhasePresence& cond = m.cond[0];

                if (Details::PhaseUsed::water(pu)) {
                    cond.setFreeWater();
                }

                if (Details::PhaseUsed::oil(pu)) {
                    cond.setFreeOil();

                    if (Details::PhaseUsed::gas(pu)) {
                        const double rsmax = attr.Rmax(io);
                        const double rs =
                            (0.0 < std::abs(in[io]))
                            ? in[ig] / in[io]
                            : (0.0 < std::abs(in[ig])) ? rsmax : 0.0;

                        if (rsmax < rs) {
                            cond.setFreeGas();
                        }

                        m.rs(0) = std::min(rs, rsmax);
                    }
                }

                if (Details::PhaseUsed::gas(pu)) {
                    if (! Details::PhaseUsed::oil(pu)) {
                        // Oil *NOT* active -- not really supported.
                        cond.setFreeGas();
                    }

                    if (Details::PhaseUsed::oil(pu)) {
                        const double rvmax = attr.Rmax(ig);
                        const double rv =
                            (0.0 < std::abs(in[ig]))
                            ? (in[io] / in[ig])
                            : (0.0 < std::abs(in[io])) ? rvmax : 0.0;

                        m.rv(0) = std::min(rv, rvmax);
                    }
                }

                return m;
            }

            typename Property::ADB
            constant(const typename Property::V& x) const
            {
                return Property::ADB::constant(x);
            }

            typename Property::ADB
            constant(const double x) const
            {
                typename Property::V y(1);
                y << x;

                return this->constant(y);
            }

            typename Property::Cells
            getRegCell(const RegionId r) const
            {
                typename Property::Cells c(1, this->attr_.cell(r));

                return c;
            }
        };
    } // namespace RateConverter
} // namespace Opm

#endif  /* OPM_RATECONVERTER_HPP_HEADER_INCLUDED */