PVT classes: implement temperature dependent viscosity

this is implemented such that if the simulation is unchanged,
viscosities are assumed to be not temperature dependent. (only if
pvtObject->set{Oil,Wat}visctTables() is called, temperature dependence
is considered.)

opm/core/props/pvt/PvtConstCompr.hpp

@@ -23,6 +23,8 @@
 #include <opm/core/props/pvt/PvtInterface.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
+
 #include <vector>
 #include <algorithm>
 
@@ -98,7 +100,9 @@ namespace Opm
               ref_B_(1, 1.0),
               comp_(1, 0.0),
               viscosity_(1, visc),
-              visc_comp_(1, 0.0)
+              visc_comp_(1, 0.0),
+              oilvisctTables_(0),
+              watvisctTables_(0)
         {
         }
 
@@ -119,6 +123,13 @@ namespace Opm
                 int tableIdx = getTableIndex_(pvtRegionIdx, i);
                 double x = -visc_comp_[tableIdx]*(p[i] - ref_press_[tableIdx]);
                 output_mu[i] = viscosity_[tableIdx]/(1.0 + x + 0.5*x*x);
+
+                if (oilvisctTables_ != 0 || watvisctTables_ != 0) {
+                    // TODO: temperature dependence
+                    OPM_THROW(std::logic_error,
+                              "temperature dependent viscosity as a function of z "
+                              "is not yet implemented!");
+                }
             }
         }
 
@@ -146,7 +157,7 @@ namespace Opm
         virtual void mu(const int n,
                         const int* pvtRegionIdx,
                         const double* p,
-                        const double* /*T*/,
+                        const double* T,
                         const double* /*r*/,
                         const PhasePresence* /*cond*/,
                         double* output_mu,
@@ -161,6 +172,46 @@ namespace Opm
                 double d = (1.0 + x + 0.5*x*x);
                 output_mu[i] = viscosity_[tableIdx]/d;
                 output_dmudp[i] = (viscosity_[tableIdx]/(d*d))*(1+x) * visc_comp_[tableIdx];
+
+                if (oilvisctTables_ != 0) {
+                    // this object handles _either_ oil _or_ water
+                    assert(watvisctTables_ == 0);
+
+                    // temperature dependence of the oil phase
+                    DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx);
+                    double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                    double muRef = -visc_comp_[tableIdx]*(pRef - ref_press_[tableIdx]);
+
+                    double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]);
+                    double alpha = muOilvisct/muRef;
+
+                    output_mu[i] *= alpha;
+                    output_dmudp[i] *= alpha;
+                    output_dmudr[i] *= alpha;
+
+                    // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                    // probably requires a healthy portion of if-spaghetti
+                }
+
+                if (watvisctTables_ != 0) {
+                    // this object handles _either_ oil _or_ water
+                    assert(oilvisctTables_ == 0);
+
+                    // temperature dependence of the water phase
+                    DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx);
+                    double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                    double muRef = -visc_comp_[tableIdx]*(pRef - ref_press_[tableIdx]);
+
+                    double muWatvisct = (*watvisctTables_)[tableIdx].evaluate("Viscosity", T[i]);
+                    double alpha = muWatvisct/muRef;
+
+                    output_mu[i] *= alpha;
+                    output_dmudp[i] *= alpha;
+                    output_dmudr[i] *= alpha;
+
+                    // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                    // probably requires a healthy portion of if-spaghetti
+                }
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
         }
@@ -287,6 +338,22 @@ namespace Opm
             std::fill(output_dRdp, output_dRdp + n, 0.0);
         }
 
+        /// set the tables which specify the temperature dependence of the oil viscosity
+        void setOilvisctTables(const std::vector<Opm::OilvisctTable>& oilvisctTables,
+                               DeckKeywordConstPtr viscrefKeyword)
+        {
+            oilvisctTables_ = &oilvisctTables;
+            viscrefKeyword_ = viscrefKeyword;
+        }
+
+        /// set the tables which specify the temperature dependence of the water viscosity
+        void setWatvisctTables(const std::vector<Opm::WatvisctTable>& watvisctTables,
+                               DeckKeywordConstPtr viscrefKeyword)
+        {
+            watvisctTables_ = &watvisctTables;
+            viscrefKeyword_ = viscrefKeyword;
+        }
+
     private:
         int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
         {
@@ -302,6 +369,10 @@ namespace Opm
         std::vector<double> comp_;
         std::vector<double> viscosity_;
         std::vector<double> visc_comp_;
+
+        const std::vector<Opm::OilvisctTable>* oilvisctTables_;
+        const std::vector<Opm::WatvisctTable>* watvisctTables_;
+        DeckKeywordConstPtr viscrefKeyword_;
     };
 
 }

opm/core/props/pvt/PvtDead.cpp

@@ -123,13 +123,20 @@ namespace Opm
             double tempInvB = b_[regionIdx](p[i]);
             double tempInvBmu = inverseBmu_[regionIdx](p[i]);
             output_mu[i] = tempInvB / tempInvBmu;
+
+            if (oilvisctTables_ != 0) {
+                // TODO: temperature dependence
+                OPM_THROW(std::logic_error,
+                          "temperature dependent viscosity as a function of z "
+                          "is not yet implemented!");
+            }
         }
     }
 
     void PvtDead::mu(const int n,
                      const int* pvtTableIdx,
                      const double* p,
-                     const double* /*T*/,
+                     const double* T,
                                const double* /*r*/,
                                double* output_mu,
                                double* output_dmudp,
@@ -143,6 +150,24 @@ namespace Opm
                 output_mu[i] = tempInvB / tempInvBmu;
                 output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
                                  - tempInvB * inverseBmu_[regionIdx].derivative(p[i])) / (tempInvBmu * tempInvBmu);
+
+                if (oilvisctTables_ != 0) {
+                    // temperature dependence of the oil phase
+                    DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx);
+                    double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                    double muRef = b_[regionIdx](pRef)/inverseBmu_[regionIdx](pRef);
+
+                    double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]);
+                    double alpha = muOilvisct/muRef;
+
+                    output_mu[i] *= alpha;
+                    output_dmudp[i] *= alpha;
+                    output_dmudr[i] *= alpha;
+
+                    // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                    // probably requires a healthy portion of if-spaghetti
+                }
+
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
 
@@ -151,7 +176,7 @@ namespace Opm
     void PvtDead::mu(const int n,
                      const int* pvtTableIdx,
                      const double* p,
-                     const double* /*T*/,
+                     const double* T,
                                const double* /*r*/,
                                const PhasePresence* /*cond*/,
                                double* output_mu,
@@ -167,6 +192,23 @@ namespace Opm
                 output_dmudp[i] = (tempInvBmu * b_[regionIdx].derivative(p[i])
                                  - tempInvB * inverseBmu_[regionIdx].derivative(p[i]))
                                  / (tempInvBmu * tempInvBmu);
+
+                if (oilvisctTables_ != 0) {
+                    // temperature dependence of the oil phase
+                    DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx);
+                    double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                    double muRef = b_[regionIdx](pRef)/inverseBmu_[regionIdx](pRef);
+
+                    double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]);
+                    double alpha = muOilvisct/muRef;
+
+                    output_mu[i] *= alpha;
+                    output_dmudp[i] *= alpha;
+                    output_dmudr[i] *= alpha;
+
+                    // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                    // probably requires a healthy portion of if-spaghetti
+                }
             }
             std::fill(output_dmudr, output_dmudr + n, 0.0);
 

opm/core/props/pvt/PvtDead.hpp

@@ -44,7 +44,11 @@ namespace Opm
     class PvtDead : public PvtInterface
     {
     public:
-        PvtDead() {};
+        PvtDead()
+        {
+            // by default, specify no temperature dependence of the PVT properties
+            oilvisctTables_ = 0;
+        }
 
         void initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables);
         void initFromGas(const std::vector<Opm::PvdgTable>& pvdgTables);
@@ -150,6 +154,15 @@ namespace Opm
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const;
+
+        /// set the tables which specify the temperature dependence of the oil viscosity
+        void setOilvisctTables(const std::vector<Opm::OilvisctTable>& oilvisctTables,
+                               DeckKeywordConstPtr viscrefKeyword)
+        {
+            oilvisctTables_ = &oilvisctTables;
+            viscrefKeyword_ = viscrefKeyword;
+        }
+
     private:
         int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
         {
@@ -163,6 +176,9 @@ namespace Opm
         std::vector<NonuniformTableLinear<double> > b_;
         std::vector<NonuniformTableLinear<double> > viscosity_;
         std::vector<NonuniformTableLinear<double> > inverseBmu_;
+
+        const std::vector<Opm::OilvisctTable>* oilvisctTables_;
+        DeckKeywordConstPtr viscrefKeyword_;
     };
 }
 

opm/core/props/pvt/PvtDeadSpline.cpp

@@ -36,7 +36,10 @@ namespace Opm
     //-------------------------------------------------------------------------
 
     PvtDeadSpline::PvtDeadSpline()
-    {}
+    {
+        // by default, specify no temperature dependence of the PVT properties
+        oilvisctTables_ = 0;
+    }
 
     void PvtDeadSpline::initFromOil(const std::vector<Opm::PvdoTable>& pvdoTables,
                                     int numSamples)
@@ -120,7 +123,7 @@ namespace Opm
     void PvtDeadSpline::mu(const int n,
                            const int* pvtTableIdx,
                            const double* p,
-                           const double* /*T*/,
+                           const double* T,
                            const double* /*r*/,
                            double* output_mu,
                            double* output_dmudp,
@@ -131,6 +134,23 @@ namespace Opm
             int regionIdx = getTableIndex_(pvtTableIdx, i);
             output_mu[i] = viscosity_[regionIdx](p[i]);
             output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
+
+            if (oilvisctTables_ != 0) {
+                // temperature dependence of the oil phase
+                DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx);
+                double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                double muRef = viscosity_[regionIdx](pRef);
+
+                double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]);
+                double alpha = muOilvisct/muRef;
+
+                output_mu[i] *= alpha;
+                output_dmudp[i] *= alpha;
+                output_dmudr[i] *= alpha;
+
+                // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                // probably requires a healthy portion of if-spaghetti
+            }
         }
         std::fill(output_dmudr, output_dmudr + n, 0.0);
     }
@@ -138,7 +158,7 @@ namespace Opm
     void PvtDeadSpline::mu(const int n,
                            const int* pvtTableIdx,
                            const double* p,
-                           const double* /*T*/,
+                           const double* T,
                            const double* /*r*/,
                            const PhasePresence* /*cond*/,
                            double* output_mu,
@@ -151,6 +171,23 @@ namespace Opm
             int regionIdx = getTableIndex_(pvtTableIdx, i);
             output_mu[i] = viscosity_[regionIdx](p[i]);
             output_dmudp[i] = viscosity_[regionIdx].derivative(p[i]);
+
+            if (oilvisctTables_ != 0) {
+                // temperature dependence of the oil phase
+                DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(regionIdx);
+                double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                double muRef = viscosity_[regionIdx](pRef);
+
+                double muOilvisct = (*oilvisctTables_)[regionIdx].evaluate("Viscosity", T[i]);
+                double alpha = muOilvisct/muRef;
+
+                output_mu[i] *= alpha;
+                output_dmudp[i] *= alpha;
+                output_dmudr[i] *= alpha;
+
+                // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                // probably requires a healthy portion of if-spaghetti
+            }
         }
         std::fill(output_dmudr, output_dmudr + n, 0.0);
     }

opm/core/props/pvt/PvtDeadSpline.hpp

@@ -148,6 +148,15 @@ namespace Opm
                           const double* z,
                           double* output_R,
                           double* output_dRdp) const;
+
+        /// set the tables which specify the temperature dependence of the oil viscosity
+        void setOilvisctTables(const std::vector<Opm::OilvisctTable>& oilvisctTables,
+                               DeckKeywordConstPtr viscrefKeyword)
+        {
+            oilvisctTables_ = &oilvisctTables;
+            viscrefKeyword_ = viscrefKeyword;
+        }
+
     private:
         int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
         {
@@ -160,6 +169,9 @@ namespace Opm
         // table per PVT region.
         std::vector<UniformTableLinear<double> > b_;
         std::vector<UniformTableLinear<double> > viscosity_;
+
+        const std::vector<Opm::OilvisctTable>* oilvisctTables_;
+        DeckKeywordConstPtr viscrefKeyword_;
     };
 
 }

opm/core/props/pvt/PvtLiveGas.cpp

@@ -114,6 +114,9 @@ namespace Opm
             double inverseBMu = miscible_gas(p[i], z + num_phases_*i, getTableIndex_(pvtRegionIdx, i), 3, false);
 
             output_mu[i] = inverseB / inverseBMu;
+
+            // temperature dependence: Since E100 does not implement temperature
+            // dependence of gas viscosity, we skip it here as well...
         }
     }
 
@@ -162,6 +165,8 @@ namespace Opm
             output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
                               / (inverseBMu * inverseBMu);
 
+            // temperature dependence: Since E100 does not implement temperature
+            // dependence of gas viscosity, we skip it here as well...
         }
 
     }

opm/core/props/pvt/PvtLiveOil.cpp

@@ -36,6 +36,9 @@ namespace Opm
     //-------------------------------------------------------------------------
     PvtLiveOil::PvtLiveOil(const std::vector<Opm::PvtoTable>& pvtoTables)
     {
+        // by default, specify no temperature dependence of the PVT properties
+        oilvisctTables_ = 0;
+
         int numTables = pvtoTables.size();
         saturated_oil_table_.resize(numTables);
         undersat_oil_tables_.resize(numTables);
@@ -136,7 +139,7 @@ namespace Opm
     void PvtLiveOil::mu(const int n,
                         const int* pvtTableIdx,
                         const double* p,
-                        const double* /*T*/,
+                        const double* T,
                               const double* z,
                               double* output_mu) const
     {
@@ -149,13 +152,20 @@ namespace Opm
 
             output_mu[i] = inverseB / inverseBMu; 
         }
+
+        if (oilvisctTables_ != 0) {
+            // TODO: temperature dependence
+            OPM_THROW(std::logic_error,
+                      "temperature dependent viscosity as a function of z "
+                      "is not yet implemented!");
+        }
     }
 
     /// Viscosity and its p and r derivatives as a function of p, T and r.
     void PvtLiveOil::mu(const int n,
                         const int* pvtTableIdx,
                         const double* p,
-                        const double* /*T*/,
+                        const double* T,
                                const double* r,
                                double* output_mu,
                                double* output_dmudp,
@@ -183,14 +193,32 @@ namespace Opm
                     output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
                                       / (inverseBMu * inverseBMu);
 
+                    if (oilvisctTables_ != 0) {
+                        // temperature dependence
+                        DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx);
+                        double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                        double rsRef = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0);
+                        double muRef = miscible_oil(pRef, rsRef, tableIdx, 3, 0)/miscible_oil(pRef, rsRef, tableIdx, 1, 0);
+
+                        double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]);
+                        double alpha = muOilvisct/muRef;
+
+                        output_mu[i] *= alpha;
+                        output_dmudp[i] *= alpha;
+                        output_dmudr[i] *= alpha;
+
+                        // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                        // probably requires a healthy portion of if-spaghetti
+                    }
                 }
+
     }
 
     /// Viscosity and its p and r derivatives as a function of p, T and r.
     void PvtLiveOil::mu(const int n,
                         const int* pvtTableIdx,
                         const double* p,
-                        const double* /*T*/,
+                        const double* T,
                                const double* r,
                                const PhasePresence* cond,
                                double* output_mu,
@@ -220,6 +248,24 @@ namespace Opm
                     output_dmudr[i] = (inverseBMu * dinverseBdr - inverseB * dinverseBmudr)
                                       / (inverseBMu * inverseBMu);
 
+
+                    if (oilvisctTables_ != 0) {
+                        // temperature dependence
+                        DeckRecordConstPtr viscrefRecord = viscrefKeyword_->getRecord(tableIdx);
+                        double pRef = viscrefRecord->getItem("REFERENCE_PRESSURE")->getSIDouble(0);
+                        double rsRef = viscrefRecord->getItem("REFERENCE_RS")->getSIDouble(0);
+                        double muRef = miscible_oil(pRef, rsRef, tableIdx, 3, 0)/miscible_oil(pRef, rsRef, tableIdx, 1, 0);
+
+                        double muOilvisct = (*oilvisctTables_)[tableIdx].evaluate("Viscosity", T[i]);
+                        double alpha = muOilvisct/muRef;
+
+                        output_mu[i] *= alpha;
+                        output_dmudp[i] *= alpha;
+                        output_dmudr[i] *= alpha;
+
+                        // TODO (?): derivative of oil viscosity w.r.t. temperature.
+                        // probably requires a healthy portion of if-spaghetti
+                    }
                 }
     }
 
@@ -645,5 +691,4 @@ namespace Opm
                 }
             }
         }
-
 } // namespace Opm

opm/core/props/pvt/PvtLiveOil.hpp

@@ -142,6 +142,14 @@ namespace Opm
                           double* output_R,
                           double* output_dRdp) const;
 
+        /// set the tables which specify the temperature dependence of the oil viscosity
+        void setOilvisctTables(const std::vector<Opm::OilvisctTable>& oilvisctTables,
+                               DeckKeywordConstPtr viscrefKeyword)
+        {
+            oilvisctTables_ = &oilvisctTables;
+            viscrefKeyword_ = viscrefKeyword;
+        }
+
     private:
         int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
         {
@@ -179,6 +187,9 @@ namespace Opm
         // store one table per PVT region.
         std::vector<std::vector<std::vector<double> > > saturated_oil_table_;
         std::vector<std::vector<std::vector<std::vector<double> > > > undersat_oil_tables_;
+
+        const std::vector<Opm::OilvisctTable>* oilvisctTables_;
+        DeckKeywordConstPtr viscrefKeyword_;
     };
 
 }