EquilibrationHelpers.hpp

Go to the documentation of this file.

/*
  Copyright 2014 SINTEF ICT, Applied Mathematics.
  Copyright 2017 IRIS
 
  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/utility/linearInterpolation.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/RootFinders.hpp>
 
#include <opm/parser/eclipse/EclipseState/InitConfig/Equil.hpp>
 
#include <opm/material/fluidsystems/BlackOilFluidSystem.hpp>
#include <opm/material/fluidstates/SimpleModularFluidState.hpp>
#include <opm/material/fluidmatrixinteractions/EclMaterialLawManager.hpp>
 
#include <memory>
 
 
/*
---- synopsis of EquilibrationHelpers.hpp ----
 
namespace Opm
{
    namespace EQUIL {
 
        namespace Miscibility {
            class RsFunction;
            class NoMixing;
            template <class FluidSystem>
            class RsVD;
            template <class FluidSystem>
            class RsSatAtContact;
        }
 
        class EquilReg;
 
 
        template <class FluidSystem,  class MaterialLaw, class MaterialLawManager>
        struct PcEq;
 
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
        inline double satFromPc(const MaterialLawManager& materialLawManager,
                                const int phase,
                                const int cell,
                                const double target_pc,
                                const bool increasing = false)
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
        inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
                                      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 {
 
 
    typedef Opm::FluidSystems::BlackOil<double> FluidSystemSimple;
 
    // Adjust oil pressure according to gas saturation and cap pressure
    typedef Opm::SimpleModularFluidState<double,
    /*numPhases=*/3,
    /*numComponents=*/3,
    FluidSystemSimple,
    /*storePressure=*/false,
    /*storeTemperature=*/false,
    /*storeComposition=*/false,
    /*storeFugacity=*/false,
    /*storeSaturation=*/true,
    /*storeDensity=*/false,
    /*storeViscosity=*/false,
    /*storeEnthalpy=*/false> SatOnlyFluidState;
 
        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;
                }
            };
 
 
            template <class FluidSystem>
            class RsVD : public RsFunction {
            public:
                RsVD(const int pvtRegionIdx,
                     const std::vector<double>& depth,
                     const std::vector<double>& rs)
                    : pvtRegionIdx_(pvtRegionIdx)
                    , depth_(depth)
                    , rs_(rs)
                {
                }
 
                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), linearInterpolationNoExtrapolation(depth_, rs_, depth));
                    }
                }
 
            private:
                const int pvtRegionIdx_;
                std::vector<double> depth_; 
                std::vector<double> rs_;    
                double satRs(const double press, const double temp) const
                {
                    return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
                }
            };
 
 
            template <class FluidSystem>
            class RvVD : public RsFunction {
            public:
                RvVD(const int pvtRegionIdx,
                     const std::vector<double>& depth,
                     const std::vector<double>& rv)
                    : pvtRegionIdx_(pvtRegionIdx)
                    , depth_(depth)
                    , rv_(rv)
                {
                }
 
                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), linearInterpolationNoExtrapolation(depth_, rv_, depth));
                    }
                }
 
            private:
                const int pvtRegionIdx_;
                std::vector<double> depth_; 
                std::vector<double> rv_;    
                double satRv(const double press, const double temp) const
                {
                    return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);
                }
            };
 
 
            template <class FluidSystem>
            class RsSatAtContact : public RsFunction {
            public:
                RsSatAtContact(const int pvtRegionIdx, const double p_contact,  const double T_contact)
                    : pvtRegionIdx_(pvtRegionIdx)
                {
                    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 int pvtRegionIdx_;
                double rs_sat_contact_;
 
                double satRs(const double press, const double temp) const
                {
                    return FluidSystem::oilPvt().saturatedGasDissolutionFactor(pvtRegionIdx_, temp, press);
                }
            };
 
 
            template <class FluidSystem>
            class RvSatAtContact : public RsFunction {
            public:
                RvSatAtContact(const int pvtRegionIdx, const double p_contact, const double T_contact)
                    :pvtRegionIdx_(pvtRegionIdx)
                {
                    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 int pvtRegionIdx_;
                double rv_sat_contact_;
 
                double satRv(const double press, const double temp) const
                {
                    return FluidSystem::gasPvt().saturatedOilVaporizationFactor(pvtRegionIdx_, temp, press);;
                }
            };
 
        } // namespace Miscibility
 
        class EquilReg {
        public:
            EquilReg(const EquilRecord& rec,
                     std::shared_ptr<Miscibility::RsFunction> rs,
                     std::shared_ptr<Miscibility::RsFunction> rv,
                     const int pvtIdx)
                : rec_    (rec)
                , rs_     (rs)
                , rv_     (rv)
                , pvtIdx_ (pvtIdx)
            {
            }
 
            typedef Miscibility::RsFunction CalcDissolution;
 
            typedef Miscibility::RsFunction CalcEvaporation;
 
            double datum()    const { return this->rec_.datumDepth(); }
 
            double pressure() const { return this->rec_.datumDepthPressure(); }
 
            double zwoc()     const { return this->rec_.waterOilContactDepth(); }
 
            double pcow_woc() const { return this->rec_.waterOilContactCapillaryPressure(); }
 
            double zgoc()     const { return this->rec_.gasOilContactDepth(); }
 
            double pcgo_goc() const { return this->rec_.gasOilContactCapillaryPressure(); }
 
 
            const CalcDissolution&
            dissolutionCalculator() const { return *this->rs_; }
 
            const CalcEvaporation&
            evaporationCalculator() const { return *this->rv_; }
 
            int pvtIdx() const { return this->pvtIdx_; }
 
 
        private:
            EquilRecord rec_;     
            std::shared_ptr<Miscibility::RsFunction> rs_;      
            std::shared_ptr<Miscibility::RsFunction> rv_;      
            const int pvtIdx_;
        };
 
 
 
        template <class FluidSystem,  class MaterialLaw, class MaterialLawManager>
        struct PcEq
        {
            PcEq(const MaterialLawManager& materialLawManager,
                 const int phase,
                 const int cell,
                 const double target_pc)
                : materialLawManager_(materialLawManager),
                  phase_(phase),
                  cell_(cell),
                  target_pc_(target_pc)
            {
 
            }
            double operator()(double s) const
            {
                const auto& matParams = materialLawManager_.materialLawParams(cell_);
                SatOnlyFluidState fluidState;
                fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
                fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
                fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
                fluidState.setSaturation(phase_, s);
 
                double pc[FluidSystem::numPhases];
                std::fill(pc, pc + FluidSystem::numPhases, 0.0);
                MaterialLaw::capillaryPressures(pc, matParams, fluidState);
                double sign = (phase_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
                double pcPhase = pc[FluidSystem::oilPhaseIdx] + sign *  pc[phase_];
                return pcPhase - target_pc_;
            }
        private:
            const MaterialLawManager& materialLawManager_;
            const int phase_;
            const int cell_;
            const double target_pc_;
        };
 
        template <class FluidSystem, class MaterialLawManager>
        double minSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) {
            const auto& scaledDrainageInfo =
                materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
 
            // Find minimum and maximum saturations.
            switch(phase) {
            case FluidSystem::waterPhaseIdx :
            {
                 return scaledDrainageInfo.Swl;
                 break;
            }
            case FluidSystem::gasPhaseIdx :
            {
                 return scaledDrainageInfo.Sgl;
                 break;
            }
            case FluidSystem::oilPhaseIdx :
            {
                 OPM_THROW(std::runtime_error, "Min saturation not implemented for oil phase.");
                 break;
            }
            default:  OPM_THROW(std::runtime_error, "Unknown phaseIdx .");
            }
            return -1.0;
        }
 
        template <class FluidSystem, class MaterialLawManager>
        double maxSaturations(const MaterialLawManager& materialLawManager, const int phase, const int cell) {
            const auto& scaledDrainageInfo =
                materialLawManager.oilWaterScaledEpsInfoDrainage(cell);
 
            // Find minimum and maximum saturations.
            switch(phase) {
            case FluidSystem::waterPhaseIdx :
            {
                 return scaledDrainageInfo.Swu;
                 break;
            }
            case FluidSystem::gasPhaseIdx :
            {
                 return scaledDrainageInfo.Sgu;
                 break;
            }
            case FluidSystem::oilPhaseIdx :
            {
                 OPM_THROW(std::runtime_error, "Max saturation not implemented for oil phase.");
                 break;
            }
            default:  OPM_THROW(std::runtime_error, "Unknown phaseIdx .");
            }
            return -1.0;
        }
 
 
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager >
        inline double satFromPc(const MaterialLawManager& materialLawManager,
                                const int phase,
                                const int cell,
                                const double target_pc,
                                const bool increasing = false)
        {
            // Find minimum and maximum saturations.           
            const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
            const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
 
            // Create the equation f(s) = pc(s) - target_pc
            const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, 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 = 60;
                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;
            }
        }
 
 
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
        struct PcEqSum
        {
            PcEqSum(const MaterialLawManager& materialLawManager,
                    const int phase1,
                    const int phase2,
                    const int cell,
                    const double target_pc)
                : materialLawManager_(materialLawManager),
                  phase1_(phase1),
                  phase2_(phase2),
                  cell_(cell),
                  target_pc_(target_pc)
            {
            }
            double operator()(double s) const
            {
                const auto& matParams = materialLawManager_.materialLawParams(cell_);
                SatOnlyFluidState fluidState;
                fluidState.setSaturation(FluidSystem::waterPhaseIdx, 0.0);
                fluidState.setSaturation(FluidSystem::oilPhaseIdx, 0.0);
                fluidState.setSaturation(FluidSystem::gasPhaseIdx, 0.0);
                fluidState.setSaturation(phase1_, s);
                fluidState.setSaturation(phase2_, 1.0 - s);
 
                double pc[FluidSystem::numPhases];
                std::fill(pc, pc + FluidSystem::numPhases, 0.0);
 
                MaterialLaw::capillaryPressures(pc, matParams, fluidState);
                double sign1 = (phase1_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
                double pc1 = pc[FluidSystem::oilPhaseIdx] + sign1 *  pc[phase1_];
                double sign2 = (phase2_ == FluidSystem::waterPhaseIdx)? -1.0 : 1.0;
                double pc2 = pc[FluidSystem::oilPhaseIdx] + sign2 *  pc[phase2_];
                return pc1 + pc2 - target_pc_;
            }
        private:
            const MaterialLawManager& materialLawManager_;
            const int phase1_;
            const int phase2_;
            const int cell_;
            const double target_pc_;
        };
 
 
 
 
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
        inline double satFromSumOfPcs(const MaterialLawManager& materialLawManager,
                                      const int phase1,
                                      const int phase2,
                                      const int cell,
                                      const double target_pc)
        {
            // Find minimum and maximum saturations.
            const double smin = minSaturations<FluidSystem>(materialLawManager, phase1, cell);
            const double smax = maxSaturations<FluidSystem>(materialLawManager, phase1, cell);
 
            // Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
            const PcEqSum<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, 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;
            }
        }
 
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
        inline double satFromDepth(const MaterialLawManager& materialLawManager,
                                   const double cellDepth,
                                   const double contactDepth,
                                   const int phase,
                                   const int cell,
                                   const bool increasing = false)
        {
            const double s0 = increasing ? maxSaturations<FluidSystem>(materialLawManager, phase, cell) : minSaturations<FluidSystem>(materialLawManager, phase, cell);
            const double s1 = increasing ? minSaturations<FluidSystem>(materialLawManager, phase, cell) : maxSaturations<FluidSystem>(materialLawManager, phase, cell);
 
            if (cellDepth < contactDepth){
                return s0;
            } else {
                return s1;
            }
 
        }
 
        template <class FluidSystem, class MaterialLaw, class MaterialLawManager>
        inline bool isConstPc(const MaterialLawManager& materialLawManager,
                              const int                          phase,
                              const int                          cell)
        {
            // Create the equation f(s) = pc(s);
            const PcEq<FluidSystem, MaterialLaw, MaterialLawManager> f(materialLawManager, phase, cell, 0);
            const double f0 = f(minSaturations<FluidSystem>(materialLawManager, phase, cell));
            const double f1 = f(maxSaturations<FluidSystem>(materialLawManager, phase, cell));
            return std::abs(f0 - f1) < std::numeric_limits<double>::epsilon();
        }
 
    } // namespace Equil
} // namespace Opm
 
 
#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED