From b1de0bf3382cd39152338ecb602b442d8c4c6d1d Mon Sep 17 00:00:00 2001 From: Andreas Lauser Date: Mon, 1 Dec 2014 18:15:20 +0100 Subject: [PATCH] unit system: add temperature this requires the possibility of specifying an offset for the SI conversion because Eclipse in its eternal wisdom chooses to specify temperatures using degrees Celsius and degrees Fahrenheit instead of using Kelvin an Rankine... --- opm/parser/eclipse/Deck/DeckDoubleItem.cpp | 13 +++----- opm/parser/eclipse/Deck/DeckFloatItem.cpp | 13 +++----- .../eclipse/Units/ConversionFactors.hpp | 17 ++++++++++ opm/parser/eclipse/Units/Dimension.cpp | 32 ++++++++++++++++--- opm/parser/eclipse/Units/Dimension.hpp | 12 +++++-- opm/parser/eclipse/Units/UnitSystem.cpp | 16 ++++++++-- opm/parser/eclipse/Units/UnitSystem.hpp | 2 +- opm/parser/eclipse/Units/tests/UnitTests.cpp | 10 ++++++ 8 files changed, 88 insertions(+), 27 deletions(-) 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/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")); }