diff --git a/opm/parser/eclipse/CMakeLists.txt b/opm/parser/eclipse/CMakeLists.txt index 476015fe7..a75e0d7a1 100644 --- a/opm/parser/eclipse/CMakeLists.txt +++ b/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 diff --git a/opm/parser/eclipse/Deck/DeckDoubleItem.cpp b/opm/parser/eclipse/Deck/DeckDoubleItem.cpp index 9445fe0bb..7ca21a9f6 100644 --- a/opm/parser/eclipse/Deck/DeckDoubleItem.cpp +++ b/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(),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_data[index] * SIfactor; - } + + for (size_t index=0; index < m_data.size(); index++) { + size_t dimIndex = (index % m_dimensions.size()); + m_SIdata[index] = m_dimensions[dimIndex]->convertRawToSi(m_data[index]); } } else throw std::invalid_argument("No dimension has been set for item:" + name() + " can not ask for SI data"); diff --git a/opm/parser/eclipse/Deck/DeckFloatItem.cpp b/opm/parser/eclipse/Deck/DeckFloatItem.cpp index d68a0f7c3..517e145e9 100644 --- a/opm/parser/eclipse/Deck/DeckFloatItem.cpp +++ b/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(),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_data[index] * SIfactor; - } + + for (size_t index=0; index < m_data.size(); index++) { + size_t dimIndex = (index % m_dimensions.size()); + m_SIdata[index] = m_dimensions[dimIndex]->convertRawToSi(m_data[index]); } } else throw std::invalid_argument("No dimension has been set for item:" + name() + " can not ask for SI data"); diff --git a/opm/parser/eclipse/EclipseState/EclipseState.cpp b/opm/parser/eclipse/EclipseState/EclipseState.cpp index 2fc856b2b..01a618cfb 100644 --- a/opm/parser/eclipse/EclipseState/EclipseState.cpp +++ b/opm/parser/eclipse/EclipseState/EclipseState.cpp @@ -200,6 +200,10 @@ namespace Opm { return m_rvvdTables; } + const std::vector& EclipseState::getRtempvdTables() const { + return m_rtempvdTables; + } + const std::vector& 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::quiet_NaN(); const auto eptLookup = std::make_shared>(*deck, *this); + const auto tempLookup = std::make_shared>(*deck, *this); const auto distributeTopLayer = std::make_shared(*this); const auto initPORV = std::make_shared(*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" ), diff --git a/opm/parser/eclipse/EclipseState/EclipseState.hpp b/opm/parser/eclipse/EclipseState/EclipseState.hpp index a6586ae4c..9ec740d0c 100644 --- a/opm/parser/eclipse/EclipseState/EclipseState.hpp +++ b/opm/parser/eclipse/EclipseState/EclipseState.hpp @@ -47,6 +47,7 @@ #include #include #include +#include #include #include @@ -104,6 +105,7 @@ namespace Opm { const std::vector& getRocktabTables() const; const std::vector& getRsvdTables() const; const std::vector& getRvvdTables() const; + const std::vector& getRtempvdTables() const; const std::vector& getSgofTables() const; const std::vector& getSwofTables() const; @@ -219,6 +221,7 @@ namespace Opm { std::vector m_rocktabTables; std::vector m_rsvdTables; std::vector m_rvvdTables; + std::vector m_rtempvdTables; std::vector m_sgofTables; std::vector m_swofTables; diff --git a/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp b/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp index a05b476da..50ea60464 100644 --- a/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp +++ b/opm/parser/eclipse/EclipseState/Grid/GridPropertyInitializers.hpp @@ -21,6 +21,7 @@ #include #include +#include #include #include @@ -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 GridPropertyTemperatureLookupInitializer + : public GridPropertyBaseInitializer +{ +public: + GridPropertyTemperatureLookupInitializer(const Deck& deck, const EclipseState& eclipseState) + : m_deck(deck) + , m_eclipseState(eclipseState) + { } + + void apply(std::vector& 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::quiet_NaN(); + std::fill(values.begin(), values.end(), nan); + return; + } + + auto eclipseGrid = m_eclipseState.getEclipseGrid(); + const std::vector& rtempvdTables = m_eclipseState.getRtempvdTables(); + const std::vector& 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 diff --git a/opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp b/opm/parser/eclipse/EclipseState/Tables/RtempvdTable.hpp new file mode 100644 index 000000000..a5db388bc --- /dev/null +++ b/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 . + */ +#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{ + "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 &getDepthColumn() const + { return ParentType::getColumn(0); } + + const std::vector &getTemperatureColumn() const + { return ParentType::getColumn(1); } + }; +} + +#endif diff --git a/opm/parser/eclipse/Units/ConversionFactors.hpp b/opm/parser/eclipse/Units/ConversionFactors.hpp index 6e1d94c9c..55d6471c4 100644 --- a/opm/parser/eclipse/Units/ConversionFactors.hpp +++ b/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; diff --git a/opm/parser/eclipse/Units/Dimension.cpp b/opm/parser/eclipse/Units/Dimension.cpp index fe325741c..0e54799e5 100644 --- a/opm/parser/eclipse/Units/Dimension.cpp +++ b/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; diff --git a/opm/parser/eclipse/Units/Dimension.hpp b/opm/parser/eclipse/Units/Dimension.hpp index 382dc3c13..acc93094f 100644 --- a/opm/parser/eclipse/Units/Dimension.hpp +++ b/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; }; } diff --git a/opm/parser/eclipse/Units/UnitSystem.cpp b/opm/parser/eclipse/Units/UnitSystem.cpp index 87b8d6dc4..6ad85b233 100644 --- a/opm/parser/eclipse/Units/UnitSystem.cpp +++ b/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) { - std::shared_ptr dim( new Dimension(dimension , SI_factor) ); + void UnitSystem::addDimension(const std::string& dimension , double SIfactor, double SIoffset) { + std::shared_ptr 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 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::newComposite( dimension , SIfactor )); @@ -106,6 +113,9 @@ namespace Opm { boost::split(parts , dimension , boost::is_any_of("/")); std::shared_ptr dividend = parseFactor( parts[0] ); std::shared_ptr 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::newComposite( dimension , dividend->getSIScaling() / divisor->getSIScaling() )); } else { @@ -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); diff --git a/opm/parser/eclipse/Units/UnitSystem.hpp b/opm/parser/eclipse/Units/UnitSystem.hpp index 1908d4b2d..d71abaf3e 100644 --- a/opm/parser/eclipse/Units/UnitSystem.hpp +++ b/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 dimension); std::shared_ptr getNewDimension(const std::string& dimension); std::shared_ptr getDimension(const std::string& dimension) const; diff --git a/opm/parser/eclipse/Units/tests/UnitTests.cpp b/opm/parser/eclipse/Units/tests/UnitTests.cpp index aeb691650..51ec0ca28 100644 --- a/opm/parser/eclipse/Units/tests/UnitTests.cpp +++ b/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 length = system.getDimension("Length"); std::shared_ptr time = system.getDimension("Time"); + std::shared_ptr 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 s BOOST_CHECK( system->hasDimension("Time")); BOOST_CHECK( system->hasDimension("Permeability")); BOOST_CHECK( system->hasDimension("Pressure")); + BOOST_CHECK( system->hasDimension("Temperature")); } diff --git a/opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD b/opm/parser/share/keywords/000_Eclipse100/R/RTEMPVD new file mode 100644 index 000000000..b98089dcd --- /dev/null +++ b/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"] + } + ] +} + + diff --git a/opm/parser/share/keywords/000_Eclipse100/T/TEMP b/opm/parser/share/keywords/000_Eclipse100/T/TEMP new file mode 100644 index 000000000..4da2ae3c6 --- /dev/null +++ b/opm/parser/share/keywords/000_Eclipse100/T/TEMP @@ -0,0 +1 @@ +{"name" : "TEMP", "sections" : ["RUNSPEC"]} diff --git a/opm/parser/share/keywords/001_Eclipse300/T/TEMPI b/opm/parser/share/keywords/001_Eclipse300/T/TEMPI new file mode 100644 index 000000000..55717c1f3 --- /dev/null +++ b/opm/parser/share/keywords/001_Eclipse300/T/TEMPI @@ -0,0 +1,8 @@ +{ + "name": "TEMPI", + "sections": ["SOLUTION"], + "data": { + "value_type": "DOUBLE", + "dimension" : "Temperature" + } +} diff --git a/opm/parser/share/keywords/001_Eclipse300/T/TEMPVD b/opm/parser/share/keywords/001_Eclipse300/T/TEMPVD new file mode 100644 index 000000000..24dc66ddc --- /dev/null +++ b/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"] + } + ] +} + +