Merge pull request #370 from andlaus/basic_temperature_support

Basic temperature support

opm/parser/eclipse/CMakeLists.txt

@@ -184,6 +184,7 @@ EclipseState/Tables/ImptvdTable.hpp
 EclipseState/Tables/RsvdTable.hpp
 EclipseState/Tables/PvtoInnerTable.hpp
 EclipseState/Tables/RvvdTable.hpp
+EclipseState/Tables/RtempvdTable.hpp
 EclipseState/Tables/PvtgOuterTable.hpp
 #
 Utility/WconinjeWrapper.hpp

opm/parser/eclipse/Deck/DeckDoubleItem.cpp

@@ -46,15 +46,10 @@ namespace Opm {
                 return;
             }
             m_SIdata.resize( m_data.size() );
-            if (m_dimensions.size() == 1) {
+
-                double SIfactor = m_dimensions[0]->getSIScaling();
-                std::transform( m_data.begin() , m_data.end() , m_SIdata.begin() , std::bind1st(std::multiplies<double>(),SIfactor));
-            } else {
             for (size_t index=0; index < m_data.size(); index++) {
                 size_t dimIndex = (index % m_dimensions.size());
-                    double SIfactor = m_dimensions[dimIndex]->getSIScaling();
+                m_SIdata[index] = m_dimensions[dimIndex]->convertRawToSi(m_data[index]);
-                    m_SIdata[index] = m_data[index] * SIfactor;
-                }
             }
         } else
             throw std::invalid_argument("No dimension has been set for item:" + name() + " can not ask for SI data");

opm/parser/eclipse/Deck/DeckFloatItem.cpp

@@ -59,15 +59,10 @@ namespace Opm {
                 return;
             }
             m_SIdata.resize( m_data.size() );
-            if (m_dimensions.size() == 1) {
+
-                float SIfactor = m_dimensions[0]->getSIScaling();
-                std::transform( m_data.begin() , m_data.end() , m_SIdata.begin() , std::bind1st(std::multiplies<float>(),SIfactor));
-            } else {
             for (size_t index=0; index < m_data.size(); index++) {
                 size_t dimIndex = (index % m_dimensions.size());
-                    float SIfactor = m_dimensions[dimIndex]->getSIScaling();
+                m_SIdata[index] = m_dimensions[dimIndex]->convertRawToSi(m_data[index]);
-                    m_SIdata[index] = m_data[index] * SIfactor;
-                }
             }
         } else
             throw std::invalid_argument("No dimension has been set for item:" + name() + " can not ask for SI data");

opm/parser/eclipse/EclipseState/EclipseState.cpp

@@ -200,6 +200,10 @@ namespace Opm {
         return m_rvvdTables;
     }
 
+    const std::vector<RtempvdTable>& EclipseState::getRtempvdTables() const {
+        return m_rtempvdTables;
+    }
+
     const std::vector<SgofTable>& EclipseState::getSgofTables() const {
         return m_sgofTables;
     }
@@ -244,6 +248,17 @@ namespace Opm {
         //depends on the presence of the RKTRMDIR keyword...
         initRocktabTables(deck, parserLog);
 
+        // the temperature vs depth table. the problem here is that
+        // the TEMPVD (E300) and RTEMPVD (E300 + E100) keywords are
+        // synonymous, but we want to provide only a single cannonical
+        // API here, so we jump through some small hoops...
+        if (deck->hasKeyword("TEMPVD") && deck->hasKeyword("RTEMPVD"))
+            throw std::invalid_argument("The TEMPVD and RTEMPVD tables are mutually exclusive!");
+        else if (deck->hasKeyword("TEMPVD"))
+            initSimpleTables(deck, parserLog, "TEMPVD", m_rtempvdTables);
+        else if (deck->hasKeyword("RTEMPVD"))
+            initSimpleTables(deck, parserLog, "RTEMPVD", m_rtempvdTables);
+
         initFullTables(deck, parserLog, "PVTG", m_pvtgTables);
         initFullTables(deck, parserLog, "PVTO", m_pvtoTables);
    }
@@ -522,6 +537,7 @@ namespace Opm {
 
         double nan = std::numeric_limits<double>::quiet_NaN();
         const auto eptLookup = std::make_shared<GridPropertyEndpointTableLookupInitializer<>>(*deck, *this);
+        const auto tempLookup = std::make_shared<GridPropertyTemperatureLookupInitializer<>>(*deck, *this);
         const auto distributeTopLayer = std::make_shared<const GridPropertyPostProcessor::DistributeTopLayer>(*this);
         const auto initPORV = std::make_shared<GridPropertyPostProcessor::InitPORV>(*this);
 
@@ -663,6 +679,9 @@ namespace Opm {
             SupportedDoubleKeywordInfo( "ISWCRZ"   , eptLookup, "1" ),
             SupportedDoubleKeywordInfo( "ISWCRZ-"  , eptLookup, "1" ),
 
+            // cell temperature (E300 only, but makes a lot of sense for E100, too)
+            SupportedDoubleKeywordInfo( "TEMPI"    , tempLookup, "Temperature" ),
+
             // porosity
             SupportedDoubleKeywordInfo( "PORO"  , nan, distributeTopLayer , "1" ),
 

opm/parser/eclipse/EclipseState/EclipseState.hpp

@@ -47,6 +47,7 @@
 #include <opm/parser/eclipse/EclipseState/Tables/RocktabTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RsvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/RvvdTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
 
@@ -104,6 +105,7 @@ namespace Opm {
         const std::vector<RocktabTable>& getRocktabTables() const;
         const std::vector<RsvdTable>& getRsvdTables() const;
         const std::vector<RvvdTable>& getRvvdTables() const;
+        const std::vector<RtempvdTable>& getRtempvdTables() const;
         const std::vector<SgofTable>& getSgofTables() const;
         const std::vector<SwofTable>& getSwofTables() const;
 
@@ -219,6 +221,7 @@ namespace Opm {
         std::vector<RocktabTable> m_rocktabTables;
         std::vector<RsvdTable> m_rsvdTables;
         std::vector<RvvdTable> m_rvvdTables;
+        std::vector<RtempvdTable> m_rtempvdTables;
         std::vector<SgofTable> m_sgofTables;
         std::vector<SwofTable> m_swofTables;
 

opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp

@@ -21,6 +21,7 @@
 
 #include <opm/parser/eclipse/EclipseState/Tables/SwofTable.hpp>
 #include <opm/parser/eclipse/EclipseState/Tables/SgofTable.hpp>
+#include <opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp>
 
 #include <vector>
 #include <string>
@@ -349,6 +350,51 @@ private:
     const Deck& m_deck;
     const EclipseState& m_eclipseState;
 };
+
+// initialize the TEMPI grid property using the temperature vs depth
+// table (stemming from the TEMPVD or the RTEMPVD keyword)
+template <class EclipseState=Opm::EclipseState,
+          class Deck=Opm::Deck>
+class GridPropertyTemperatureLookupInitializer
+    : public GridPropertyBaseInitializer<double>
+{
+public:
+    GridPropertyTemperatureLookupInitializer(const Deck& deck, const EclipseState& eclipseState)
+        : m_deck(deck)
+        , m_eclipseState(eclipseState)
+    { }
+
+    void apply(std::vector<double>& values,
+               const std::string& propertyName) const
+    {
+        if (propertyName != "TEMPI")
+            throw std::logic_error("The TemperatureLookupInitializer can only be used for the initial temperature!");
+
+        if (!m_deck.hasKeyword("EQLNUM")) {
+            // if values are defaulted in the TEMPI keyword, but no
+            // EQLNUM is specified, you will get NaNs...
+            double nan = std::numeric_limits<double>::quiet_NaN();
+            std::fill(values.begin(), values.end(), nan);
+            return;
+        }
+
+        auto eclipseGrid = m_eclipseState.getEclipseGrid();
+        const std::vector<RtempvdTable>& rtempvdTables = m_eclipseState.getRtempvdTables();
+        const std::vector<int>& eqlNum = m_eclipseState.getIntGridProperty("EQLNUM")->getData();
+
+        for (int cellIdx = 0; cellIdx < eqlNum.size(); ++ cellIdx) {
+            int cellEquilNum = eqlNum[cellIdx];
+            const RtempvdTable& rtempvdTable = rtempvdTables[cellEquilNum];
+            double cellDepth = std::get<2>(eclipseGrid->getCellCenter(cellIdx));
+            values[cellIdx] = rtempvdTable.evaluate("Temperature", cellDepth);
+        }
+    }
+
+private:
+    const Deck& m_deck;
+    const EclipseState& m_eclipseState;
+};
+
 }
 
 #endif

opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp

@@ -0,0 +1,68 @@
+/*
+  Copyright (C) 2014 by Andreas Lauser
+
+  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_PARSER_RTEMPVD_TABLE_HPP
+#define	OPM_PARSER_RTEMPVD_TABLE_HPP
+
+#include "SingleRecordTable.hpp"
+
+namespace Opm {
+    // forward declaration
+    class EclipseState;
+
+    class RtempvdTable : protected SingleRecordTable {
+        typedef SingleRecordTable ParentType;
+
+        friend class EclipseState;
+        RtempvdTable() = default;
+
+        /*!
+         * \brief Read the RTEMPVD keyword and provide some convenience
+         *        methods for it.
+         */
+        void init(Opm::DeckKeywordConstPtr keyword, int recordIdx)
+        {
+            ParentType::init(keyword,
+                             std::vector<std::string>{
+                                 "Depth",
+                                 "Temperature"
+                             },
+                             recordIdx,
+                             /*firstEntityOffset=*/0);
+
+            ParentType::checkNonDefaultable("Depth");
+            ParentType::checkMonotonic("Depth", /*isAscending=*/true);
+
+            ParentType::checkNonDefaultable("Temperature");
+        }
+
+    public:
+        using ParentType::numTables;
+        using ParentType::numRows;
+        using ParentType::numColumns;
+        using ParentType::evaluate;
+
+        const std::vector<double> &getDepthColumn() const
+        { return ParentType::getColumn(0); }
+
+        const std::vector<double> &getTemperatureColumn() const
+        { return ParentType::getColumn(1); }
+    };
+}
+
+#endif

opm/parser/eclipse/Units/ConversionFactors.hpp

@@ -143,6 +143,19 @@ namespace Opm {
             const double psia   = lbf / square(inch);
             /// @}
 
+            /// \name Temperature. This one is more complicated
+            /// because the unit systems used by Eclipse (i.e. degrees
+            /// Celsius and degrees Fahrenheit require to add or
+            /// subtract an offset for the conversion between from/to
+            /// Kelvin
+            /// @{
+            const double degCelsius = 1.0; // scaling factor °C -> K
+            const double degCelsiusOffset = 273.15; // offset for the °C -> K conversion
+
+            const double degFahrenheit = 5.0/9; // scaling factor °F -> K
+            const double degFahrenheitOffset = 255.37; // offset for the °C -> K conversion
+            /// @}
+
             /// \name Viscosity
             /// @{
             const double Pas   = Pascal * second; // == 1
@@ -179,6 +192,8 @@ namespace Opm {
         using namespace details::prefix;
         using namespace details::unit;
         const double Pressure             = barsa;
+        const double Temperature          = degCelsius;
+        const double TemperatureOffset    = degCelsiusOffset;
         const double Length               = meter;
         const double Time                 = day;
         const double Mass                 = kilogram;
@@ -198,6 +213,8 @@ namespace Opm {
         using namespace details::prefix;
         using namespace details::unit;
         const double Pressure             = psia;
+        const double Temperature          = degFahrenheit;
+        const double TemperatureOffset    = degFahrenheitOffset;
         const double Length               = feet;
         const double Time                 = day;
         const double Mass                 = pound;

opm/parser/eclipse/Units/Dimension.cpp

@@ -29,7 +29,7 @@ namespace Opm {
         
     }
 
-    Dimension::Dimension(const std::string& name, double SIfactor) 
+    Dimension::Dimension(const std::string& name, double SIfactor, double SIoffset)
     {
         for (auto iter = name.begin(); iter != name.end(); ++iter) {
             if (!isalpha(*iter) && (*iter) != '1')
@@ -37,9 +37,9 @@ namespace Opm {
         }
         m_name = name;
         m_SIfactor = SIfactor;
+        m_SIoffset = SIoffset;
     }
 
-    
     double Dimension::getSIScaling() const {
         if (!std::isfinite(m_SIfactor))
             throw std::logic_error("The DeckItem contains a field with a context dependent unit. "
@@ -47,15 +47,39 @@ namespace Opm {
         return m_SIfactor;
     }
 
+    double Dimension::getSIOffset() const {
+        return m_SIoffset;
+    }
+
+    double Dimension::convertRawToSi(double rawValue) const {
+        if (!std::isfinite(m_SIfactor))
+            throw std::logic_error("The DeckItem contains a field with a context dependent unit. "
+                                   "Use getRawDoubleData() and convert the returned value manually!");
+
+        return rawValue*m_SIfactor + m_SIoffset;
+    }
+
+    double Dimension::convertSiToRaw(double siValue) const {
+        if (!std::isfinite(m_SIfactor))
+            throw std::logic_error("The DeckItem contains a field with a context dependent unit. "
+                                   "Use getRawDoubleData() and convert the returned value manually!");
+
+        return (siValue - m_SIoffset)/m_SIfactor;
+    }
+
     const std::string& Dimension::getName() const {
         return m_name;
     }
 
+    // only dimensions with zero offset are compositable...
+    bool Dimension::isCompositable() const
+    { return m_SIoffset == 0.0; }
 
-    Dimension * Dimension::newComposite(const std::string& dim , double SIfactor) {
+    Dimension * Dimension::newComposite(const std::string& dim , double SIfactor, double SIoffset) {
         Dimension * dimension = new Dimension();
         dimension->m_name = dim;
         dimension->m_SIfactor = SIfactor;
+        dimension->m_SIoffset = SIoffset;
 
         return dimension;
     }
@@ -64,7 +88,7 @@ namespace Opm {
     bool Dimension::equal(const Dimension& other) const {
         if (m_name != other.m_name)
             return false;
-        if (m_SIfactor == other.m_SIfactor)
+        if (m_SIfactor == other.m_SIfactor && m_SIoffset == other.m_SIoffset)
             return true;
         if (std::isnan(m_SIfactor) && std::isnan(other.m_SIfactor))
             return true;

opm/parser/eclipse/Units/Dimension.hpp

@@ -26,16 +26,24 @@ namespace Opm {
 
     class Dimension {
     public:
-        Dimension(const std::string& name, double SI_factor);
+        Dimension(const std::string& name, double SIfactor, double SIoffset = 0.0);
+
         double getSIScaling() const;
+        double getSIOffset() const;
+
+        double convertRawToSi(double rawValue) const;
+        double convertSiToRaw(double siValue) const;
+
         bool equal(const Dimension& other) const;
         const std::string& getName() const;
-        static Dimension * newComposite(const std::string& dim , double SIfactor);
+        bool isCompositable() const;
+        static Dimension * newComposite(const std::string& dim, double SIfactor, double SIoffset = 0.0);
 
     private:
         Dimension();
         std::string m_name;
         double m_SIfactor;
+        double m_SIoffset;
     };
 }
 

opm/parser/eclipse/Units/UnitSystem.cpp

@@ -66,8 +66,8 @@ namespace Opm {
         m_dimensions.insert( std::make_pair(dimension->getName() , dimension));
     }
 
-    void UnitSystem::addDimension(const std::string& dimension , double SI_factor) {
+    void UnitSystem::addDimension(const std::string& dimension , double SIfactor, double SIoffset) {
-        std::shared_ptr<const Dimension> dim( new Dimension(dimension , SI_factor) );
+        std::shared_ptr<const Dimension> dim( new Dimension(dimension , SIfactor, SIoffset) );
         addDimension(dim);
     }
 
@@ -82,6 +82,13 @@ namespace Opm {
         double SIfactor = 1.0;
         for (auto iter = dimensionList.begin(); iter != dimensionList.end(); ++iter) {
             std::shared_ptr<const Dimension> dim = getDimension( *iter );
+
+            // all constituing dimension must be compositable. The
+            // only exception is if there is the "composite" dimension
+            // consists of exactly a single atomic dimension...
+            if (dimensionList.size() > 1 && !dim->isCompositable())
+                throw std::invalid_argument("Composite dimensions currently cannot require a conversion offset");
+
             SIfactor *= dim->getSIScaling();
         }
         return std::shared_ptr<Dimension>(Dimension::newComposite( dimension , SIfactor ));
@@ -107,6 +114,9 @@ namespace Opm {
             std::shared_ptr<const Dimension> dividend = parseFactor( parts[0] );
             std::shared_ptr<const Dimension> divisor = parseFactor( parts[1] );
 
+            if (dividend->getSIOffset() != 0.0 || divisor->getSIOffset() != 0.0)
+                throw std::invalid_argument("Composite dimensions cannot currently require a conversion offset");
+        
             return std::shared_ptr<Dimension>( Dimension::newComposite( dimension , dividend->getSIScaling() / divisor->getSIScaling() ));
         } else {
             return parseFactor( dimension );
@@ -139,6 +149,7 @@ namespace Opm {
         
         system->addDimension("1"         , 1.0);
         system->addDimension("Pressure"  , Metric::Pressure );
+        system->addDimension("Temperature", Metric::Temperature, Metric::TemperatureOffset);
         system->addDimension("Length"    , Metric::Length);
         system->addDimension("Time"      , Metric::Time );
         system->addDimension("Mass"         , Metric::Mass );
@@ -162,6 +173,7 @@ namespace Opm {
         
         system->addDimension("1"    , 1.0);
         system->addDimension("Pressure", Field::Pressure );
+        system->addDimension("Temperature", Field::Temperature, Field::TemperatureOffset);
         system->addDimension("Length", Field::Length);
         system->addDimension("Time" , Field::Time);
         system->addDimension("Mass", Field::Mass);

opm/parser/eclipse/Units/UnitSystem.hpp

@@ -33,7 +33,7 @@ namespace Opm {
         UnitSystem(const std::string& unitSystem);
         const std::string& getName() const;
 
-        void addDimension(const std::string& dimension , double SI_factor);
+        void addDimension(const std::string& dimension, double SIfactor, double SIoffset = 0.0);
         void addDimension(std::shared_ptr<const Dimension> dimension);
         std::shared_ptr<const Dimension> getNewDimension(const std::string& dimension);
         std::shared_ptr<const Dimension> getDimension(const std::string& dimension) const;

opm/parser/eclipse/Units/tests/UnitTests.cpp

@@ -41,6 +41,11 @@ BOOST_AUTO_TEST_CASE(makeComposite) {
     BOOST_CHECK_EQUAL(100 , composite->getSIScaling());
 }
 
+BOOST_AUTO_TEST_CASE(MakeCompositeInvalid) {
+    // conversion to SI temperatures requires an offset, but such
+    // composite units are (currently?) invalid
+    BOOST_CHECK_THROW(Dimension::newComposite("Length*Temperature", 100), std::logic_error);
+}
 
 BOOST_AUTO_TEST_CASE(CreateDimensionInvalidNameThrows) {
     BOOST_CHECK_THROW(Dimension(" " , 1) , std::invalid_argument);
@@ -89,11 +94,15 @@ BOOST_AUTO_TEST_CASE(UnitSystemAddDimensions) {
     UnitSystem system("Metric");
     system.addDimension("Length" , 1 );
     system.addDimension("Time" , 86400 );
+    system.addDimension("Temperature", 1.0, 273.15);
 
     std::shared_ptr<const Dimension> length = system.getDimension("Length");
     std::shared_ptr<const Dimension> time = system.getDimension("Time");
+    std::shared_ptr<const Dimension> temperature = system.getDimension("Temperature");
     BOOST_CHECK_EQUAL(1     , length->getSIScaling());
     BOOST_CHECK_EQUAL(86400 , time->getSIScaling());
+    BOOST_CHECK_EQUAL(1.0   , temperature->getSIScaling());
+    BOOST_CHECK_EQUAL(273.15, temperature->getSIOffset());
 
     system.addDimension("Length" , 0.3048);
     length = system.getDimension("Length");
@@ -123,6 +132,7 @@ static void checkSystemHasRequiredDimensions(std::shared_ptr<const UnitSystem> s
     BOOST_CHECK( system->hasDimension("Time"));
     BOOST_CHECK( system->hasDimension("Permeability"));
     BOOST_CHECK( system->hasDimension("Pressure"));
+    BOOST_CHECK( system->hasDimension("Temperature"));
 }
 
 

opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD

@@ -0,0 +1,18 @@
+{
+    "name" : "RTEMPVD",
+    "sections" : [ "PROPS", "SOLUTION" ],
+    "size" : {
+       "keyword" : "EQLDIMS" ,
+       "item" : "NTEQUL"
+    },
+    "items" : [
+      {
+         "name":"DATA",
+         "value_type":"DOUBLE",
+          "size_type" : "ALL",
+          "dimension" : ["Length","Temperature"]
+      }
+    ]
+}
+
+

opm/parser/share/keywords/000_Eclipse100/T/TEMP

@@ -0,0 +1 @@
+{"name" : "TEMP", "sections" : ["RUNSPEC"]}

opm/parser/share/keywords/001_Eclipse300/T/TEMPI

@@ -0,0 +1,8 @@
+{
+ "name": "TEMPI",
+ "sections": ["SOLUTION"],
+ "data": {
+   "value_type": "DOUBLE",
+   "dimension" : "Temperature"
+ }
+}

opm/parser/share/keywords/001_Eclipse300/T/TEMPVD

@@ -0,0 +1,18 @@
+{
+    "name" : "TEMPVD",
+    "sections" : [ "PROPS", "SOLUTION" ],
+    "size" : {
+       "keyword" : "EQLDIMS" ,
+       "item" : "NTEQUL"
+    },
+    "items" : [
+      {
+         "name":"DATA",
+         "value_type":"DOUBLE",
+          "size_type" : "ALL",
+          "dimension" : ["Length","Temperature"]
+      }
+    ]
+}
+
+