diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 55adafcb2..f221a492b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -22,10 +22,12 @@ list (APPEND MAIN_SOURCE_FILES opm/common/data/SimulationDataContainer.cpp + opm/common/util/numeric/cmp.cpp ) list (APPEND TEST_SOURCE_FILES tests/test_SimulationDataContainer.cpp + tests/test_cmp.cpp ) list (APPEND TEST_DATA_FILES @@ -42,6 +44,7 @@ list (APPEND PROGRAM_SOURCE_FILES list( APPEND PUBLIC_HEADER_FILES opm/common/data/SimulationDataContainer.hpp + opm/common/util/numeric/cmp.hpp opm/common/ErrorMacros.hpp opm/common/Exceptions.hpp opm/common/utility/platform_dependent/disable_warnings.h diff --git a/opm/common/data/SimulationDataContainer.cpp b/opm/common/data/SimulationDataContainer.cpp index 524da544f..9cabf925b 100644 --- a/opm/common/data/SimulationDataContainer.cpp +++ b/opm/common/data/SimulationDataContainer.cpp @@ -17,7 +17,7 @@ along with OPM. If not, see . */ - +#include #include namespace Opm { @@ -52,15 +52,26 @@ namespace Opm { std::vector& SimulationDataContainer::getCellData( const std::string& name ) { auto iter = m_cell_data.find( name ); - if (iter == m_cell_data.end()) + if (iter == m_cell_data.end()) { throw std::invalid_argument("The cell data with name: " + name + " does not exist"); - else + } else return iter->second; } + + const std::vector& SimulationDataContainer::getCellData( const std::string& name ) const { + auto iter = m_cell_data.find( name ); + if (iter == m_cell_data.end()) { + throw std::invalid_argument("The cell data with name: " + name + " does not exist"); + } else + return iter->second; + } + + void SimulationDataContainer::registerCellData( const std::string& name , size_t components , double initialValue) { - if (!hasCellData( name )) + if (!hasCellData( name )) { m_cell_data.insert( std::pair>( name , std::vector(components * m_num_cells , initialValue ))); + } } @@ -71,18 +82,66 @@ namespace Opm { std::vector& SimulationDataContainer::getFaceData( const std::string& name ) { auto iter = m_face_data.find( name ); - if (iter == m_face_data.end()) - throw std::invalid_argument("The Face data with name: " + name + " does not exist"); - else + if (iter == m_face_data.end()) { + throw std::invalid_argument("The face data with name: " + name + " does not exist"); + } else return iter->second; } - void SimulationDataContainer::registerFaceData( const std::string& name , size_t components , double initialValue) { - if (!hasFaceData( name )) - m_face_data.insert( std::pair>( name , std::vector(components * m_num_faces , initialValue ))); + const std::vector& SimulationDataContainer::getFaceData( const std::string& name ) const { + auto iter = m_face_data.find( name ); + if (iter == m_face_data.end()) { + throw std::invalid_argument("The Face data with name: " + name + " does not exist"); + } else + return iter->second; } + void SimulationDataContainer::registerFaceData( const std::string& name , size_t components , double initialValue) { + if (!hasFaceData( name )) { + m_face_data.insert( std::pair>( name , std::vector(components * m_num_faces , initialValue ))); + } + } + + bool SimulationDataContainer::equal( const SimulationDataContainer& other ) const { + if ((m_num_cells != other.m_num_cells) || + (m_num_phases != other.m_num_phases) || + (m_num_faces != other.m_num_faces)) + return false; + + if ((m_face_data.size() != other.m_face_data.size()) || + (m_cell_data.size() != other.m_cell_data.size())) + return false; + + for (const auto& cell_data : m_cell_data) { + const auto key = cell_data.first; + const auto data = cell_data.second; + + if (other.hasCellData( key )) { + const auto& other_data = other.getCellData( key ); + if (!cmp::double_vector_equal( data , other_data )) + return false; + } else + return false; + } + + for (const auto& face_data : m_face_data) { + const auto key = face_data.first; + const auto data = face_data.second; + + if (other.hasFaceData( key )) { + const auto& other_data = other.getFaceData( key ); + if (!cmp::double_vector_equal( data , other_data )) + return false; + } else + return false; + } + + return true; + } + + + /* This is very deprecated. */ void SimulationDataContainer::addDefaultFields() { registerCellData("PRESSURE" , 1 , 0.0); diff --git a/opm/common/data/SimulationDataContainer.hpp b/opm/common/data/SimulationDataContainer.hpp index d2d98f9e8..51f1fb396 100644 --- a/opm/common/data/SimulationDataContainer.hpp +++ b/opm/common/data/SimulationDataContainer.hpp @@ -27,6 +27,21 @@ namespace Opm { + + /// The SimulationDataContainer is a simple container to manage + /// simulation data. The container is instantiated with information + /// of how many cells, faces and phases are present in the + /// reservoirmodel. You can then add data to the container by using the + /// + /// registerCellData() + /// registerFaceData() + /// + /// functions. The container owns and manages the data, but + /// mutable references are returned with the getCellData() and + /// getFaceData() methods, and the content will typically be + /// modified by external scope. + + class SimulationDataContainer { public: SimulationDataContainer(size_t num_cells, size_t num_faces , size_t num_phases); @@ -36,12 +51,19 @@ namespace Opm { size_t numCells() const; bool hasCellData( const std::string& name ) const; + + /// Will register a data vector of size numCells() * + /// components. void registerCellData( const std::string& name , size_t components , double initialValue = 0.0 ); std::vector& getCellData( const std::string& name ); + const std::vector& getCellData( const std::string& name ) const; bool hasFaceData( const std::string& name ) const; void registerFaceData( const std::string& name , size_t components , double initialValue = 0.0 ); std::vector& getFaceData( const std::string& name ); + const std::vector& getFaceData( const std::string& name ) const; + + bool equal(const SimulationDataContainer& other) const; /* Old deprecated */ diff --git a/opm/common/util/numeric/cmp.cpp b/opm/common/util/numeric/cmp.cpp new file mode 100644 index 000000000..7d77aa24b --- /dev/null +++ b/opm/common/util/numeric/cmp.cpp @@ -0,0 +1,91 @@ +/* + Copyright 2016 Statoil ASA. + + 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 . +*/ + +#include +#include + +#include +#include +#include + +#include + + +namespace Opm { + namespace cmp { + + bool double_equal(double value1 , double value2 , double abs_eps , double rel_eps) { + bool equal = true; + double diff = fabs(value1 - value2); + if (diff > abs_eps) { + double scale = std::max(fabs(value1), fabs(value2)); + + if (diff > scale * rel_eps) { + equal = false; + } + } + return equal; + } + + + bool double_equal(double value1 , double value2) { + return double_equal( value1 , value2 , default_abs_epsilon , default_rel_epsilon ); + } + + /*****************************************************************/ + + bool double_vector_equal(const std::vector& v1, const std::vector& v2, double abs_eps, double rel_eps) { + if (v1.size() != v2.size()) { + return false; + } + + for (size_t i = 0; i < v1.size(); i++) { + if (!double_equal( v1[i] , v2[i])) + return false; + } + + return true; + } + + bool double_vector_equal(const std::vector& v1, const std::vector& v2) { + return double_vector_equal( v1 , v2 , default_abs_epsilon , default_rel_epsilon); + } + + + bool double_ptr_equal(const double * p1, const double *p2, size_t num_elements, double abs_eps , double rel_eps) { + if (memcmp(p1 , p2 , num_elements * sizeof * p1) == 0) + return true; + else { + size_t index; + for (index = 0; index < num_elements; index++) { + if (!double_equal( p1[index] , p2[index] , abs_eps , rel_eps)) { + return false; + } + } + } + return true; + } + + + bool double_ptr_equal(const double * p1, const double *p2, size_t num_elements) { + return double_ptr_equal( p1 , p2 ,num_elements , default_abs_epsilon , default_rel_epsilon ); + } + + } +} diff --git a/opm/common/util/numeric/cmp.hpp b/opm/common/util/numeric/cmp.hpp new file mode 100644 index 000000000..127f0e4be --- /dev/null +++ b/opm/common/util/numeric/cmp.hpp @@ -0,0 +1,67 @@ +/* + Copyright 2016 Statoil ASA. + + 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 COMMON_UTIL_NUMERIC_CMP +#define COMMON_UTIL_NUMERIC_CMP + +#include +#include + +namespace Opm { + + /// In the namespace cmp are implemented functions for + /// approximate comparison of double values based on absolute + /// and relative difference. There are three functions: + /// + /// double_equal() : Compare two double values. + /// + /// double_ptr_equal(): This compares all the element in the + /// two double * pointers. + /// + /// double_vector_equal(): This compares all the elements in + /// two std::vector instances. + /// + /// For both double_vector_equal() and double_ptr_equal() the + /// actual comparison is based on the double_equal() + /// function. All functions exist as two overloads, one which + /// takes explicit input values for the absolute and relative + /// epsilon, and one which uses default values. + /// + /// For more details of floating point comparison please consult + /// this reference: + /// + /// https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/ + namespace cmp { + + const double default_abs_epsilon = 1e-8; + const double default_rel_epsilon = 1e-5; + + + bool double_equal(double value1, double value2, double abs_eps , double rel_eps); + bool double_equal(double value1, double value2); + + bool double_vector_equal(const std::vector& v1, const std::vector& v2, double abs_eps, double rel_eps); + bool double_vector_equal(const std::vector& v1, const std::vector& v2); + + bool double_ptr_equal(const double* p1, const double* p2, size_t num_elements, double abs_eps, double rel_eps); + bool double_ptr_equal(const double* p1, const double* p2, size_t num_elements); + } +} + +#endif diff --git a/tests/test_SimulationDataContainer.cpp b/tests/test_SimulationDataContainer.cpp index 73baab131..3aaf803b1 100644 --- a/tests/test_SimulationDataContainer.cpp +++ b/tests/test_SimulationDataContainer.cpp @@ -113,3 +113,63 @@ BOOST_AUTO_TEST_CASE(TestRegisterCellData) { } +BOOST_AUTO_TEST_CASE(Test_Equal) { + { + SimulationDataContainer container1(100 , 10 , 2); + SimulationDataContainer container2(100 , 10 , 2); + BOOST_CHECK( container1.equal( container2 )); + } + + { + SimulationDataContainer container1(100 , 10 , 2); + SimulationDataContainer container2(100 , 10 , 1); + BOOST_CHECK( !container1.equal( container2 )); + } + + { + SimulationDataContainer container1(100 , 10 , 2); + SimulationDataContainer container2(100 , 10 , 2); + + container1.registerCellData( "FIELDX" , 1 , 123 ); + BOOST_CHECK( !container1.equal( container2 )); + container2.registerCellData( "FIELDX" , 1 , 123 ); + BOOST_CHECK( container1.equal( container2 )); + + container1.registerFaceData( "FACEX" , 1 , 123 ); + BOOST_CHECK( !container1.equal( container2 )); + container2.registerFaceData( "FACEX" , 1 , 123 ); + BOOST_CHECK( container1.equal( container2 )); + } + + { + SimulationDataContainer container1(100 , 10 , 2); + SimulationDataContainer container2(100 , 10 , 2); + + container1.registerCellData( "FIELD1" , 1 , 123 ); + container2.registerCellData( "FIELD2" , 1 , 123 ); + BOOST_CHECK( !container1.equal( container2 )); + } + + { + SimulationDataContainer container1(100 , 10 , 2); + SimulationDataContainer container2(100 , 10 , 2); + + container1.registerFaceData( "FIELD1" , 1 , 123 ); + container2.registerFaceData( "FIELD2" , 1 , 123 ); + BOOST_CHECK( !container1.equal( container2 )); + } + + { + SimulationDataContainer container1(100 , 10 , 2); + SimulationDataContainer container2(100 , 10 , 2); + + container1.registerFaceData( "FIELD1" , 1 , 123 ); + container2.registerFaceData( "FIELD1" , 1 , 123 ); + BOOST_CHECK( container1.equal( container2 )); + + std::vector& f = container1.getFaceData( "FIELD1" ); + f[0] *= 1.1; + BOOST_CHECK( !container1.equal( container2 )); + } +} + diff --git a/tests/test_cmp.cpp b/tests/test_cmp.cpp new file mode 100644 index 000000000..e93647437 --- /dev/null +++ b/tests/test_cmp.cpp @@ -0,0 +1,109 @@ +/* + Copyright 2016 Statoil ASA. + + 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 . + */ + + +#define BOOST_TEST_DYN_LINK +#define BOOST_TEST_MODULE FLOAT_CMP_TESTS +#include + +#include + +#include + +using namespace Opm; + +/** + Ahhh - the joys of comparing floating point numbers .... + + http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm +*/ + +BOOST_AUTO_TEST_CASE(TestSCalarcmp) { + const double abs_epsilon = cmp::default_abs_epsilon; + const double rel_epsilon = cmp::default_rel_epsilon; + + BOOST_CHECK( cmp::double_equal(1,1)); + BOOST_CHECK_EQUAL( false , cmp::double_equal(1,0)); + BOOST_CHECK_EQUAL( false , cmp::double_equal(0,1)); + BOOST_CHECK_EQUAL( false , cmp::double_equal(-1,1)); + + + + double v1,v2; + /* Should be equal: */ + { + v1 = 0.0; + v2 = 0.0; + BOOST_CHECK( cmp::double_equal( v1 , v2)); + + v1 = 1e-12; + v2 = v1 + 0.5*abs_epsilon; + BOOST_CHECK( cmp::double_equal( v1 , v2)); + + v1 = 7.0; + v2 = 7.0; + BOOST_CHECK( cmp::double_equal( v1 , v2)); + + v1 = -7.0; + v2 = -7.0; + BOOST_CHECK( cmp::double_equal( v1 , v2)); + + v1 = 0; + v2 = 0.5 * abs_epsilon; + BOOST_CHECK( cmp::double_equal( v1 , v2)); + + + v1 = 1e7; + v2 = 1e7 + 2*abs_epsilon; + BOOST_CHECK( cmp::double_equal( v1 , v2 )); + + v1 = 1e7*(1 - abs_epsilon); + v2 = 1e7*(1 + rel_epsilon); + BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + + v1 = 1e7*(1 + abs_epsilon); + v2 = 1e7*(1 + rel_epsilon); + BOOST_CHECK( cmp::double_equal( v1 , v2 )); + } + + /* Should be different: */ + { + v1 = 0; + v2 = 1.5 * abs_epsilon; + BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + + v1 = 1e-8; + v2 = v1 + 1.5*abs_epsilon; + BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + + v1 = 1; + v2 = v1*(1 + 2*rel_epsilon + abs_epsilon); + BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + + v1 = 10; + v2 = v1*(1 + 2*rel_epsilon + abs_epsilon); + BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + + v1 = 1e7; + v2 = 1e7*(1 + 2*rel_epsilon + abs_epsilon); + BOOST_CHECK( !cmp::double_equal( v1 , v2 )); + } +} + +