Merge pull request #3513 from goncalvesmachadoc/gasvisct

Implement GASVISCT for black-oil [Bug fix]

opm/input/eclipse/EclipseState/Tables/GasvisctTable.hpp

@@ -24,15 +24,14 @@
 
 namespace Opm {
 
-    class Deck;
     class DeckItem;
 
     class GasvisctTable : public SimpleTable {
     public:
-        GasvisctTable( const Deck& deck, const DeckItem& deckItem );
+        GasvisctTable( const DeckItem& item, const int tableID  );
 
         const TableColumn& getTemperatureColumn() const;
-        const TableColumn& getGasViscosityColumn(size_t compIdx) const;
+        const TableColumn& getGasViscosityColumn() const;
     };
 }
 

opm/input/eclipse/EclipseState/Tables/TableManager.hpp

@@ -289,7 +289,6 @@ namespace Opm {
         void initRTempTables(const Deck& deck);
         void initDims(const Deck& deck);
         void initRocktabTables(const Deck& deck);
-        void initGasvisctTables(const Deck& deck);
 
         void initPlymaxTables(const Deck& deck);
         void initPlyrockTables(const Deck& deck);

opm/material/fluidsystems/blackoilpvt/GasPvtThermal.hpp

@@ -57,11 +57,14 @@ public:
         enableThermalDensity_ = false;
         enableJouleThomson_ = false;
         enableThermalViscosity_ = false;
+        enableInternalEnergy_ = false;
         isothermalPvt_ = nullptr;
     }
 
     GasPvtThermal(IsothermalPvt* isothermalPvt,
                   const std::vector<TabulatedOneDFunction>& gasvisctCurves,
+                  const std::vector<Scalar>& viscrefPress,
+                  const std::vector<Scalar>& viscRef,
                   const std::vector<Scalar>& gasdentRefTemp,
                   const std::vector<Scalar>& gasdentCT1,
                   const std::vector<Scalar>& gasdentCT2,
@@ -74,6 +77,8 @@ public:
                   bool enableInternalEnergy)
         : isothermalPvt_(isothermalPvt)
         , gasvisctCurves_(gasvisctCurves)
+        , viscrefPress_(viscrefPress)
+        , viscRef_(viscRef)
         , gasdentRefTemp_(gasdentRefTemp)
         , gasdentCT1_(gasdentCT1)
         , gasdentCT2_(gasdentCT2)
@@ -105,6 +110,8 @@ public:
     void setNumRegions(size_t numRegions)
     {
         gasvisctCurves_.resize(numRegions);
+        viscrefPress_.resize(numRegions);
+        viscRef_.resize(numRegions);
         internalEnergyCurves_.resize(numRegions);
         gasdentRefTemp_.resize(numRegions);
         gasdentCT1_.resize(numRegions);
@@ -120,9 +127,6 @@ public:
     void initEnd()
     { }
 
-    size_t numRegions() const
-    { return gasvisctCurves_.size(); }
-
     /*!
      * \brief Returns true iff the density of the gas phase is temperature dependent.
      */
@@ -141,6 +145,10 @@ public:
     bool enableThermalViscosity() const
     { return enableThermalViscosity_; }
 
+    size_t numRegions() const
+    { return viscrefPress_.size(); }
+ 
+
     /*!
      * \brief Returns the specific internal energy [J/kg] of gas given a set of parameters.
      */
@@ -216,12 +224,14 @@ public:
                          const Evaluation& Rv,
                          const Evaluation& Rvw) const
     {
+        
+        const auto& isothermalMu = isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw);
         if (!enableThermalViscosity())
-            return isothermalPvt_->viscosity(regionIdx, temperature, pressure, Rv, Rvw);
+            return isothermalMu;
 
-        // compute the viscosity deviation due to temperature
+        // compute the viscosity deviation due to temperature        
-        const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature);
+        const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, /*extrapolate=*/true);
-        return muGasvisct;
+        return muGasvisct/viscRef_[regionIdx]*isothermalMu;
     }
 
     /*!
@@ -232,12 +242,13 @@ public:
                                   const Evaluation& temperature,
                                   const Evaluation& pressure) const
     {
+        const auto& isothermalMu = isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
         if (!enableThermalViscosity())
-            return isothermalPvt_->saturatedViscosity(regionIdx, temperature, pressure);
+            return isothermalMu;
 
         // compute the viscosity deviation due to temperature
         const auto& muGasvisct = gasvisctCurves_[regionIdx].eval(temperature, true);
-        return muGasvisct;
+        return muGasvisct/viscRef_[regionIdx]*isothermalMu;
     }
 
     /*!
@@ -379,6 +390,12 @@ public:
     const std::vector<TabulatedOneDFunction>& gasvisctCurves() const
     { return gasvisctCurves_; }
 
+    const std::vector<Scalar>& viscrefPress() const
+    { return viscrefPress_; }
+
+    const std::vector<Scalar>& viscRef() const
+    { return viscRef_; }
+
     const std::vector<Scalar>& gasdentRefTemp() const
     { return gasdentRefTemp_; }
 
@@ -409,6 +426,8 @@ public:
             return false;
 
         return  this->gasvisctCurves() == data.gasvisctCurves() &&
+                this->viscrefPress() == data.viscrefPress() &&
+                this->viscRef() == data.viscRef() &&
                 this->gasdentRefTemp() == data.gasdentRefTemp() &&
                 this->gasdentCT1() == data.gasdentCT1() &&
                 this->gasdentCT2() == data.gasdentCT2() &&
@@ -428,6 +447,8 @@ public:
         else
             isothermalPvt_ = nullptr;
         gasvisctCurves_ = data.gasvisctCurves_;
+        viscrefPress_ = data.viscrefPress_;
+        viscRef_ = data.viscRef_;
         gasdentRefTemp_ = data.gasdentRefTemp_;
         gasdentCT1_ = data.gasdentCT1_;
         gasdentCT2_ = data.gasdentCT2_;
@@ -448,6 +469,8 @@ private:
     // The PVT properties needed for temperature dependence of the viscosity. We need
     // to store one value per PVT region.
     std::vector<TabulatedOneDFunction> gasvisctCurves_;
+    std::vector<Scalar> viscrefPress_;
+    std::vector<Scalar> viscRef_;
 
     std::vector<Scalar> gasdentRefTemp_;
     std::vector<Scalar> gasdentCT1_;

opm/material/fluidsystems/blackoilpvt/WaterPvtThermal.hpp

@@ -55,6 +55,7 @@ public:
     WaterPvtThermal()
     {
         enableThermalDensity_ = false;
+        enableJouleThomson_ = false;
         enableThermalViscosity_ = false;
         enableInternalEnergy_ = false;
         isothermalPvt_ = nullptr;

src/opm/input/eclipse/EclipseState/Tables/TableManager.cpp

@@ -581,6 +581,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
         initSimpleTableContainer<SpecheatTable>(deck, "SPECHEAT", m_tabdims.getNumPVTTables());
         initSimpleTableContainer<SpecrockTable>(deck, "SPECROCK", m_tabdims.getNumSatTables());
         initSimpleTableContainer<OilvisctTable>(deck, "OILVISCT", m_tabdims.getNumPVTTables());
+        initSimpleTableContainer<GasvisctTable>(deck, "GASVISCT", m_tabdims.getNumPVTTables());
         initSimpleTableContainer<WatvisctTable>(deck, "WATVISCT", m_tabdims.getNumPVTTables());
 
         initSimpleTableContainer<PlyadsTable>(deck, "PLYADS", m_tabdims.getNumSatTables());
@@ -592,7 +593,6 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
 
         initPlyrockTables(deck);
         initPlymaxTables(deck);
-        initGasvisctTables(deck);
         initRTempTables(deck);
         initRocktabTables(deck);
         initPlyshlogTables(deck);
@@ -618,34 +618,7 @@ std::optional<JFunc> make_jfunc(const Deck& deck) {
         } else if (hasRTEMPVD) {
             initSimpleTableContainer<RtempvdTable>(deck, "RTEMPVD", "RTEMPVD", m_eqldims.getNumEquilRegions());
         }
-    }
+    }   
-
-
-    void TableManager::initGasvisctTables(const Deck& deck) {
-
-        const std::string keywordName = "GASVISCT";
-        size_t numTables = m_tabdims.getNumPVTTables();
-
-        if (!deck.hasKeyword(keywordName))
-            return; // the table is not featured by the deck...
-
-        auto& container = forceGetTables(keywordName , numTables);
-
-        if (deck.count(keywordName) > 1) {
-            complainAboutAmbiguousKeyword(deck, keywordName);
-            return;
-        }
-
-        const auto& tableKeyword = deck[keywordName].back();
-        for (size_t tableIdx = 0; tableIdx < tableKeyword.size(); ++tableIdx) {
-            const auto& tableRecord = tableKeyword.getRecord( tableIdx );
-            const auto& dataItem = tableRecord.getItem( 0 );
-            if (dataItem.data_size() > 0) {
-                std::shared_ptr<GasvisctTable> table = std::make_shared<GasvisctTable>( deck , dataItem );
-                container.addTable( tableIdx , table );
-            }
-        }
-    }
 
 
     void TableManager::initPlyshlogTables(const Deck& deck) {

src/opm/input/eclipse/EclipseState/Tables/Tables.cpp

@@ -1170,46 +1170,12 @@ WatvisctTable::getWaterViscosityColumn() const
     return SimpleTable::getColumn(1);
 }
 
-GasvisctTable::GasvisctTable(const Deck& deck, const DeckItem& deckItem)
+GasvisctTable::GasvisctTable(const DeckItem& item, const int tableID)
 {
-    int numComponents = deck.get<ParserKeywords::COMPS>().back().getRecord(0).getItem(0).get<int>(0);
-
-    auto temperatureDimension = deck.getActiveUnitSystem().getDimension("Temperature");
-    auto viscosityDimension = deck.getActiveUnitSystem().getDimension("Viscosity");
 
     m_schema.addColumn(ColumnSchema("Temperature", Table::STRICTLY_INCREASING, Table::DEFAULT_NONE));
-    for (int compIdx = 0; compIdx < numComponents; ++compIdx) {
+    m_schema.addColumn(ColumnSchema("Viscosity", Table::RANDOM, Table::DEFAULT_NONE));
-        std::string columnName = "Viscosity" + std::to_string(compIdx);
+    SimpleTable::init("GASVISCT", item, tableID);
-        m_schema.addColumn(ColumnSchema(columnName, Table::INCREASING, Table::DEFAULT_NONE));
-    }
-
-    SimpleTable::addColumns();
-
-    if (deckItem.data_size() % numColumns() != 0)
-        throw std::runtime_error("Number of columns in the data file is inconsistent "
-                                 "with the expected number for keyword GASVISCT");
-
-    size_t rows = deckItem.data_size() / m_schema.size();
-    for (size_t columnIndex = 0; columnIndex < m_schema.size(); columnIndex++) {
-        auto& column = getColumn(columnIndex);
-        for (size_t rowIdx = 0; rowIdx < rows; rowIdx++) {
-            size_t deckIndex = rowIdx * m_schema.size() + columnIndex;
-
-            if (deckItem.defaultApplied(deckIndex))
-                column.addDefault("GASVISCT");
-            else {
-                double rawValue = deckItem.get<double>(deckIndex);
-                double SIValue;
-
-                if (columnIndex == 0)
-                    SIValue = temperatureDimension.convertRawToSi(rawValue);
-                else
-                    SIValue = viscosityDimension.convertRawToSi(rawValue);
-
-                column.addValue(SIValue, "GASVISCT");
-            }
-        }
-    }
 }
 
 const TableColumn&
@@ -1219,17 +1185,16 @@ GasvisctTable::getTemperatureColumn() const
 }
 
 const TableColumn&
-GasvisctTable::getGasViscosityColumn(size_t compIdx) const
+GasvisctTable::getGasViscosityColumn() const
 {
-    return SimpleTable::getColumn(1 + compIdx);
+    return SimpleTable::getColumn(1);
 }
 
 RtempvdTable::RtempvdTable(const DeckItem& item, const int tableID)
 {
     m_schema.addColumn(ColumnSchema("Depth", Table::STRICTLY_INCREASING, Table::DEFAULT_NONE));
     m_schema.addColumn(ColumnSchema("Temperature", Table::RANDOM, Table::DEFAULT_NONE));
-
+    SimpleTable::init("RTEMPVD", item, tableID);
-    SimpleTable::init("GASVISCT", item, tableID);
 }
 
 const TableColumn&

src/opm/material/fluidsystems/blackoilpvt/GasPvtThermal.cpp

@@ -61,14 +61,45 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
 
     // viscosity
     if (enableThermalViscosity_) {
+        if (tables.getViscrefTable().empty())
+            OPM_THROW(std::runtime_error, "VISCREF is required when GASVISCT is present");
+
         const auto& gasvisctTables = tables.getGasvisctTables();
-        auto gasCompIdx = tables.gas_comp_index();
+        const auto& viscrefTable = tables.getViscrefTable();
-        std::string gasvisctColumnName = "Viscosity" + std::to_string(static_cast<long long>(gasCompIdx));
+
+        if (gasvisctTables.size() != numRegions) {
+            OPM_THROW(std::runtime_error,
+                      fmt::format("Tables sizes mismatch. GASVISCT: {}, NumRegions: {}\n",
+                                  gasvisctTables.size(), numRegions));
+        }
+        if (viscrefTable.size() != numRegions) {
+            OPM_THROW(std::runtime_error,
+                      fmt::format("Tables sizes mismatch. VISCREF: {}, NumRegions: {}\n",
+                                  viscrefTable.size(), numRegions));
+        }
 
         for (unsigned regionIdx = 0; regionIdx < numRegions; ++regionIdx) {
-            const auto& T = gasvisctTables[regionIdx].getColumn("Temperature").vectorCopy();
+            const auto& TCol = gasvisctTables[regionIdx].getColumn("Temperature").vectorCopy();
-            const auto& mu = gasvisctTables[regionIdx].getColumn(gasvisctColumnName).vectorCopy();
+            const auto& muCol = gasvisctTables[regionIdx].getColumn("Viscosity").vectorCopy();
-            gasvisctCurves_[regionIdx].setXYContainers(T, mu);
+            gasvisctCurves_[regionIdx].setXYContainers(TCol, muCol);
+
+            viscrefPress_[regionIdx] = viscrefTable[regionIdx].reference_pressure;
+
+            // Temperature used to calculate the reference viscosity [K]. 
+            // The value dooes not matter as the underlying PVT object is isothermal.
+            constexpr const Scalar Tref = 273.15 + 20;
+            //TODO: For now I just assumed the default references RV and RVW = 0, 
+            // we could add these two parameters to a new item keyword VISCREF 
+            //or create a new keyword for gas. 
+            constexpr const Scalar Rvref = 0.0;
+            constexpr const Scalar Rvwref = 0.0;
+            // compute the reference viscosity using the isothermal PVT object.
+            viscRef_[regionIdx] =
+                isothermalPvt_->viscosity(regionIdx,
+                                          Tref,
+                                          viscrefPress_[regionIdx],
+                                          Rvref, 
+                                          Rvwref);
         }
     }
 

src/opm/material/fluidsystems/blackoilpvt/OilPvtThermal.cpp

@@ -52,6 +52,7 @@ initFromState(const EclipseState& eclState, const Schedule& schedule)
     const auto& tables = eclState.getTableManager();
 
     enableThermalDensity_ = tables.OilDenT().size() > 0;
+    enableJouleThomson_ = tables.OilJT().size() > 0;
     enableThermalViscosity_ = tables.hasTables("OILVISCT");
     enableInternalEnergy_ = tables.hasTables("SPECHEAT");