EquilibrationHelpers.hpp

Go to the documentation of this file.

/*
  Copyright 2014 SINTEF ICT, Applied Mathematics.

  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_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED

#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/RootFinders.hpp>

#include <memory>


/*
---- synopsis of EquilibrationHelpers.hpp ----

namespace Opm
{
    namespace Equil {

        template <class Props>
        class DensityCalculator;

        template <>
        class DensityCalculator< BlackoilPropertiesInterface >;

        namespace Miscibility {
            class RsFunction;
            class NoMixing;
            class RsVD;
            class RsSatAtContact;
        }

        struct EquilRecord;

        template <class DensCalc>
        class EquilReg;


        struct PcEq;

        inline double satFromPc(const BlackoilPropertiesInterface& props,
                                const int phase,
                                const int cell,
                                const double target_pc,
                                const bool increasing = false);
        struct PcEqSum
        inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
                                      const int phase1,
                                      const int phase2,
                                      const int cell,
                                      const double target_pc);
    } // namespace Equil
} // namespace Opm

---- end of synopsis of EquilibrationHelpers.hpp ----
*/


namespace Opm
{
    namespace Equil {


        template <class Props>
        class DensityCalculator;

        template <>
        class DensityCalculator< BlackoilPropertiesInterface > {
        public:
            DensityCalculator(const BlackoilPropertiesInterface& props,
                              const int                          c)
                : props_(props)
                , c_(1, c)
            {
            }

            std::vector<double>
            operator()(const double               p,
                       const double               T,
                       const std::vector<double>& z) const
            {
                const int np = props_.numPhases();
                std::vector<double> A(np * np, 0);

                assert (z.size() == std::vector<double>::size_type(np));

                double* dAdp = 0;
                props_.matrix(1, &p, &T, &z[0], &c_[0], &A[0], dAdp);

                std::vector<double> rho(np, 0.0);
                props_.density(1, &A[0], &c_[0], &rho[0]);

                return rho;
            }

        private:
            const BlackoilPropertiesInterface& props_;
            const std::vector<int>             c_;
        };


        namespace Miscibility {

            class RsFunction
            {
            public:
                virtual double operator()(const double depth,
                                          const double press,
                                          const double temp,
                                          const double sat = 0.0) const = 0;
            };


            class NoMixing : public RsFunction {
            public:
                double
                operator()(const double /* depth */,
                           const double /* press */,
                           const double /* temp */,
                           const double /* sat */ = 0.0) const
                {
                    return 0.0;
                }
            };


            class RsVD : public RsFunction {
            public:
                RsVD(const BlackoilPropertiesInterface& props,
                     const int cell,
                     const std::vector<double>& depth,
                     const std::vector<double>& rs)
                    : props_(props) 
                    , cell_(cell)
                    , depth_(depth)
                    , rs_(rs)
                {
                    auto pu = props_.phaseUsage();
                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
                }

                double
                operator()(const double depth,
                           const double press,
                           const double temp,
                           const double sat_gas = 0.0) const
                {
                    if (sat_gas > 0.0) {
                        return satRs(press, temp);
                    } else {
                        return std::min(satRs(press, temp), linearInterpolation(depth_, rs_, depth));
                    }
                }

            private:
                const BlackoilPropertiesInterface& props_;
                const int cell_;
                std::vector<double> depth_; 
                std::vector<double> rs_;    
                double z_[BlackoilPhases::MaxNumPhases];
                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];

                double satRs(const double press, const double temp) const
                {
                    props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
                    // Rs/Bo is in the gas row and oil column of A_.
                    // 1/Bo is in the oil row and column.
                    // Recall also that it is stored in column-major order.
                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
                    const int np = props_.numPhases();
                    return A_[np*opos + gpos] / A_[np*opos + opos];
                }
            };


            class RvVD : public RsFunction {
            public:
                RvVD(const BlackoilPropertiesInterface& props,
                     const int cell,
                     const std::vector<double>& depth,
                     const std::vector<double>& rv)
                    : props_(props) 
                    , cell_(cell)
                    , depth_(depth)
                    , rv_(rv)
                {
                    auto pu = props_.phaseUsage();
                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
                }

                double
                operator()(const double depth,
                           const double press,
                           const double temp,
                           const double sat_oil = 0.0 ) const
                {
                    if (std::abs(sat_oil) > 1e-16) {
                        return satRv(press, temp);
                    } else {
                        return std::min(satRv(press, temp), linearInterpolation(depth_, rv_, depth));
                    }
                }

            private:
                const BlackoilPropertiesInterface& props_;
                const int cell_;
                std::vector<double> depth_; 
                std::vector<double> rv_;    
                double z_[BlackoilPhases::MaxNumPhases];
                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];

                double satRv(const double press, const double temp) const
                {
                    props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
                    // Rv/Bg is in the oil row and gas column of A_.
                    // 1/Bg is in the gas row and column.
                    // Recall also that it is stored in column-major order.
                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
                    const int np = props_.numPhases();
                    return A_[np*gpos + opos] / A_[np*gpos + gpos];
                }
            };


            class RsSatAtContact : public RsFunction {
            public:
                RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact,  const double T_contact)
                    : props_(props), cell_(cell)
                {
                    auto pu = props_.phaseUsage();
                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
                    rs_sat_contact_ = satRs(p_contact, T_contact);
                }

                double
                operator()(const double /* depth */,
                           const double press,
                           const double temp,
                           const double sat_gas = 0.0) const
                {                     
                    if (sat_gas > 0.0) {
                        return satRs(press, temp);
                    } else {
                        return std::min(satRs(press, temp), rs_sat_contact_);
                    }
                }

            private:
                const BlackoilPropertiesInterface& props_;
                const int cell_;
                double z_[BlackoilPhases::MaxNumPhases];
                double rs_sat_contact_;
                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];

                double satRs(const double press, const double temp) const
                {
                    props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
                    // Rs/Bo is in the gas row and oil column of A_.
                    // 1/Bo is in the oil row and column.
                    // Recall also that it is stored in column-major order.
                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
                    const int np = props_.numPhases();
                    return A_[np*opos + gpos] / A_[np*opos + opos];
                }
            };


            class RvSatAtContact : public RsFunction {
            public:
                RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact, const double T_contact)
                    : props_(props), cell_(cell)
                {
                    auto pu = props_.phaseUsage();
                    std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
                    z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
                    z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
                    rv_sat_contact_ = satRv(p_contact, T_contact);
                }

                double
                operator()(const double /*depth*/,
                           const double press,
                           const double temp,
                           const double sat_oil = 0.0) const
                {
                    if (sat_oil > 0.0) {
                        return satRv(press, temp);
                    } else {
                        return std::min(satRv(press, temp), rv_sat_contact_);
                    }
                }

            private:
                const BlackoilPropertiesInterface& props_;
                const int cell_;
                double z_[BlackoilPhases::MaxNumPhases];
                double rv_sat_contact_;
                mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];

                double satRv(const double press, const double temp) const
                {
                    props_.matrix(1, &press, &temp, z_, &cell_, A_, 0);
                    // Rv/Bg is in the oil row and gas column of A_.
                    // 1/Bg is in the gas row and column.
                    // Recall also that it is stored in column-major order.
                    const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
                    const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
                    const int np = props_.numPhases();
                    return A_[np*gpos + opos] / A_[np*gpos + gpos];
                }
            };

        } // namespace Miscibility



        struct EquilRecord {
            struct {
                double depth;
                double press;
            } main, woc, goc;
            int live_oil_table_index;
            int wet_gas_table_index;
            int N;
        };

        template <class DensCalc>
        class EquilReg {
        public:
            EquilReg(const EquilRecord& rec,
                     const DensCalc&    density,
                     std::shared_ptr<Miscibility::RsFunction> rs,
                     std::shared_ptr<Miscibility::RsFunction> rv,
                     const PhaseUsage&  pu)
                : rec_    (rec)
                , density_(density)
                , rs_     (rs)
                , rv_     (rv)
                , pu_     (pu)
            {
            }

            typedef DensCalc CalcDensity;

            typedef Miscibility::RsFunction CalcDissolution;

            typedef Miscibility::RsFunction CalcEvaporation;

            double datum()    const { return this->rec_.main.depth; }

            double pressure() const { return this->rec_.main.press; }

            double zwoc()     const { return this->rec_.woc .depth; }

            double pcow_woc() const { return this->rec_.woc .press; }

            double zgoc()     const { return this->rec_.goc .depth; }

            double pcgo_goc() const { return this->rec_.goc .press; }

            const CalcDensity&
            densityCalculator() const { return this->density_; }

            const CalcDissolution&
            dissolutionCalculator() const { return *this->rs_; }

            const CalcEvaporation&
            evaporationCalculator() const { return *this->rv_; }

            const PhaseUsage&
            phaseUsage() const { return this->pu_; }

        private:
            EquilRecord rec_;     
            DensCalc    density_; 
            std::shared_ptr<Miscibility::RsFunction> rs_;      
            std::shared_ptr<Miscibility::RsFunction> rv_;      
            PhaseUsage  pu_;      
        };



        struct PcEq
        {
            PcEq(const BlackoilPropertiesInterface& props,
                 const int phase,
                 const int cell,
                 const double target_pc)
                : props_(props),
                  phase_(phase),
                  cell_(cell),
                  target_pc_(target_pc)
            {
                std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
                std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
            }
            double operator()(double s) const
            {
                s_[phase_] = s;
                props_.capPress(1, s_, &cell_, pc_, 0);
                return pc_[phase_] - target_pc_;
            }
        private:
            const BlackoilPropertiesInterface& props_;
            const int phase_;
            const int cell_;
            const int target_pc_;
            mutable double s_[BlackoilPhases::MaxNumPhases];
            mutable double pc_[BlackoilPhases::MaxNumPhases];
        };



        inline double satFromPc(const BlackoilPropertiesInterface& props,
                                const int phase,
                                const int cell,
                                const double target_pc,
                                const bool increasing = false)
        {
            // Find minimum and maximum saturations.
            double sminarr[BlackoilPhases::MaxNumPhases];
            double smaxarr[BlackoilPhases::MaxNumPhases];
            props.satRange(1, &cell, sminarr, smaxarr);
            const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
            const double s1 = increasing ? sminarr[phase] : smaxarr[phase];

            // Create the equation f(s) = pc(s) - target_pc
            const PcEq f(props, phase, cell, target_pc);
            const double f0 = f(s0);
            const double f1 = f(s1);

            if (f0 <= 0.0) {
                return s0;
            } else if (f1 > 0.0) {
                return s1;
            } else {
                const int max_iter = 30;
                const double tol = 1e-6;
                int iter_used = -1;
                typedef RegulaFalsi<ThrowOnError> ScalarSolver;
                const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
                return sol;
            }
        }


        struct PcEqSum
        {
            PcEqSum(const BlackoilPropertiesInterface& props,
                    const int phase1,
                    const int phase2,
                    const int cell,
                    const double target_pc)
                : props_(props),
                  phase1_(phase1),
                  phase2_(phase2),
                  cell_(cell),
                  target_pc_(target_pc)
            {
                std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
                std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
            }
            double operator()(double s) const
            {
                s_[phase1_] = s;
                s_[phase2_] = 1.0 - s;
                props_.capPress(1, s_, &cell_, pc_, 0);
                return pc_[phase1_] + pc_[phase2_] - target_pc_;
            }
        private:
            const BlackoilPropertiesInterface& props_;
            const int phase1_;
            const int phase2_;
            const int cell_;
            const int target_pc_;
            mutable double s_[BlackoilPhases::MaxNumPhases];
            mutable double pc_[BlackoilPhases::MaxNumPhases];
        };




        inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
                                      const int phase1,
                                      const int phase2,
                                      const int cell,
                                      const double target_pc)
        {
            // Find minimum and maximum saturations.
            double sminarr[BlackoilPhases::MaxNumPhases];
            double smaxarr[BlackoilPhases::MaxNumPhases];
            props.satRange(1, &cell, sminarr, smaxarr);
            const double smin = sminarr[phase1];
            const double smax = smaxarr[phase1];

            // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
            const PcEqSum f(props, phase1, phase2, cell, target_pc);
            const double f0 = f(smin);
            const double f1 = f(smax);
            if (f0 <= 0.0) {
                return smin;
            } else if (f1 > 0.0) {
                return smax;
            } else {
                const int max_iter = 30;
                const double tol = 1e-6;
                int iter_used = -1;
                typedef RegulaFalsi<ThrowOnError> ScalarSolver;
                const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
                return sol;
            }
        }

        inline double satFromDepth(const BlackoilPropertiesInterface& props,
                                   const double cellDepth,
                                   const double contactDepth,
                                   const int phase,
                                   const int cell,
                                   const bool increasing = false)
        {
            // Find minimum and maximum saturations.
            double sminarr[BlackoilPhases::MaxNumPhases];
            double smaxarr[BlackoilPhases::MaxNumPhases];
            props.satRange(1, &cell, sminarr, smaxarr);
            const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
            const double s1 = increasing ? sminarr[phase] : smaxarr[phase];

            if (cellDepth < contactDepth){
                return s0;
            } else {
                return s1;
            }

        }

        inline bool isConstPc(const BlackoilPropertiesInterface& props,
                              const int                          phase,
                              const int                          cell)
        {
            // Find minimum and maximum saturations.
            double sminarr[BlackoilPhases::MaxNumPhases];
            double smaxarr[BlackoilPhases::MaxNumPhases];
            props.satRange(1, &cell, sminarr, smaxarr);

            // Create the equation f(s) = pc(s);
            const PcEq f(props, phase, cell, 0);
            const double f0 = f(sminarr[phase]);
            const double f1 = f(smaxarr[phase]);
            return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
        }

    } // namespace Equil
} // namespace Opm


#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED