MiscibilityWater.hpp

Go to the documentation of this file.

//===========================================================================
//
// File: MiscibilityWater.hpp
//
// Created: Tue May 18 10:26:13 2010
//
// Author(s): Atgeirr F Rasmussen <atgeirr@sintef.no>
//            Bjørn Spjelkavik    <bsp@sintef.no>
//
// $Date$
//
// $Revision$
//
//===========================================================================
/*
  Copyright 2010 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 OPENRS_MISCIBILITYWATER_HEADER
#define OPENRS_MISCIBILITYWATER_HEADER

#include "MiscibilityProps.hpp"
#include <opm/common/ErrorMacros.hpp>

// Forward declaration.
class PVTW;

namespace Opm
{
    class MiscibilityWater : public MiscibilityProps
    {
    public:
        typedef std::vector<std::vector<double> > table_t;

        MiscibilityWater(DeckKeywordConstPtr pvtwKeyword)
        {
            auto pvtwRecord = pvtwKeyword->getRecord(0);
            ref_press_ = pvtwRecord->getItem("P_REF")->getSIDouble(0);
            ref_B_     = pvtwRecord->getItem("WATER_VOL_FACTOR")->getSIDouble(0);
            comp_      = pvtwRecord->getItem("WATER_COMPRESSIBILITY")->getSIDouble(0);
            viscosity_ = pvtwRecord->getItem("WATER_VISCOSITY")->getSIDouble(0);

            if (pvtwRecord->getItem("WATER_VISCOSIBILITY")->getSIDouble(0) != 0.0) {
                OPM_THROW(std::runtime_error, "MiscibilityWater does not support 'viscosibility'.");
            }
        }


        MiscibilityWater(double visc)
            : ref_press_(0.0),
              ref_B_(1.0),
              comp_(0.0),
              viscosity_(visc)
        {
        }

        // WTF?? we initialize a class for water from a keyword for oil?
        void initFromPvcdo(DeckKeywordConstPtr pvcdoKeyword)
        {
            auto pvcdoRecord = pvcdoKeyword->getRecord(0);
            ref_press_ = pvcdoRecord->getItem("P_REF")->getSIDouble(0);
            ref_B_     = pvcdoRecord->getItem("OIL_VOL_FACTOR")->getSIDouble(0);
            comp_      = pvcdoRecord->getItem("OIL_COMPRESSIBILITY")->getSIDouble(0);
            viscosity_ = pvcdoRecord->getItem("OIL_VISCOSITY")->getSIDouble(0);
            if (pvcdoRecord->getItem("OIL_VISCOSIBILITY")->getSIDouble(0) != 0.0) {
                OPM_THROW(std::runtime_error, "MiscibilityWater does not support 'viscosibility'.");
            }
        }

    virtual ~MiscibilityWater()
        {
        }

        virtual double getViscosity(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
        {
            return viscosity_;
        }

        virtual void getViscosity(const std::vector<PhaseVec>& pressures,
                                  const std::vector<CompVec>&,
                                  int,
                                  std::vector<double>& output) const
        {
            int num = pressures.size();
            output.clear();
            output.resize(num, viscosity_);
        }

        virtual double B(int /*region*/, double press, const surfvol_t& /*surfvol*/) const
        {
            if (comp_) {
                // Computing a polynomial approximation to the exponential.
                double x = comp_*(press - ref_press_);
                return ref_B_/(1.0 + x + 0.5*x*x);
            } else {
                return ref_B_;
            }
        }

        virtual void B(const std::vector<PhaseVec>& pressures,
                       const std::vector<CompVec>&,
                       int phase,
                       std::vector<double>& output) const
        {
            int num = pressures.size();
            if (comp_) {
                output.resize(num);
#pragma omp parallel for
                for (int i = 0; i < num; ++i) {
                    // Computing a polynomial approximation to the exponential.
                    double x = comp_*(pressures[i][phase] - ref_press_);
                    output[i] = ref_B_/(1.0 + x + 0.5*x*x);
                }
            } else {
                output.clear();
                output.resize(num, ref_B_);
            }
        }

    virtual double dBdp(int region, double press, const surfvol_t& surfvol) const
        {
            if (comp_) {
                return -comp_*B(region, press, surfvol);
            } else {
                return 0.0;
            }
        }

        virtual void dBdp(const std::vector<PhaseVec>& pressures,
                          const std::vector<CompVec>& surfvols,
                          int phase,
                          std::vector<double>& output_B,
                          std::vector<double>& output_dBdp) const
        {
            B(pressures, surfvols, phase, output_B);
            int num = pressures.size();
            if (comp_) {
                output_dBdp.resize(num);
#pragma omp parallel for
                for (int i = 0; i < num; ++i) {
                    output_dBdp[i] = -comp_*output_B[i];
                }
            } else {
                output_dBdp.clear();
                output_dBdp.resize(num, 0.0);
            }
        }

    virtual double R(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
        {
            return 0.0;
        }

        virtual void R(const std::vector<PhaseVec>& pressures,
                       const std::vector<CompVec>&,
                       int,
                       std::vector<double>& output) const
        {
            int num = pressures.size();
            output.clear();
            output.resize(num, 0.0);
        }

    virtual double dRdp(int /*region*/, double /*press*/, const surfvol_t& /*surfvol*/) const
        {
            return 0.0;
        }

        virtual void dRdp(const std::vector<PhaseVec>& pressures,
                          const std::vector<CompVec>&,
                          int, 
                          std::vector<double>& output_R,
                          std::vector<double>& output_dRdp) const
        {
            int num = pressures.size();
            output_R.clear();
            output_R.resize(num, 0.0);
            output_dRdp.clear();
            output_dRdp.resize(num, 0.0);
        }

    private:
        double ref_press_;
        double ref_B_;
        double comp_;
        double viscosity_;
    };

}

#endif // OPENRS_MISCIBILITYWATER_HEADER