Merge pull request #74 from joakim-hove/float-cmp

Added util/numeric/cmp.hpp

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

opm/common/data/SimulationDataContainer.cpp

@@ -17,7 +17,7 @@
   along with OPM.  If not, see <http://www.gnu.org/licenses/>.
  */
 
-
+#include <opm/common/util/numeric/cmp.hpp>
 #include <opm/common/data/SimulationDataContainer.hpp>
 
 namespace Opm {
@@ -52,15 +52,26 @@ namespace Opm {
 
     std::vector<double>& 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<double>& 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<std::string , std::vector<double>>( name , std::vector<double>(components * m_num_cells , initialValue )));
+        }
     }
 
 
@@ -71,18 +82,66 @@ namespace Opm {
 
     std::vector<double>& SimulationDataContainer::getFaceData( const std::string& name ) {
         auto iter = m_face_data.find( name );
-        if (iter == m_face_data.end())
+        if (iter == m_face_data.end()) {
-            throw std::invalid_argument("The Face data with name: " + name + " does not exist");
+            throw std::invalid_argument("The face data with name: " + name + " does not exist");
-        else
+        } else
             return iter->second;
     }
 
-    void SimulationDataContainer::registerFaceData( const std::string& name , size_t components , double initialValue) {
+    const std::vector<double>& SimulationDataContainer::getFaceData( const std::string& name ) const {
-        if (!hasFaceData( name ))
+        auto iter = m_face_data.find( name );
-            m_face_data.insert( std::pair<std::string , std::vector<double>>( name , std::vector<double>(components * m_num_faces , initialValue )));
+        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<std::string , std::vector<double>>( name , std::vector<double>(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);

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<double>& getCellData( const std::string& name );
+        const std::vector<double>& 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<double>& getFaceData( const std::string& name );
+        const std::vector<double>& getFaceData( const std::string& name ) const;
+
+        bool equal(const SimulationDataContainer& other) const;
 
 
         /* Old deprecated */

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 <http://www.gnu.org/licenses/>.
+*/
+
+#include <math.h>
+#include <string.h>
+
+#include <algorithm>
+#include <cstddef>
+#include <vector>
+
+#include <opm/common/util/numeric/cmp.hpp>
+
+
+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<double>& v1, const std::vector<double>& 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<double>& v1, const std::vector<double>& 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 );
+        }
+
+    }
+}

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 <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef COMMON_UTIL_NUMERIC_CMP
+#define COMMON_UTIL_NUMERIC_CMP
+
+#include <cstddef>
+#include <vector>
+
+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<double> 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<double>& v1, const std::vector<double>& v2, double abs_eps, double rel_eps);
+        bool double_vector_equal(const std::vector<double>& v1, const std::vector<double>& 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

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<double>& f = container1.getFaceData( "FIELD1" );
+        f[0] *= 1.1;
+        BOOST_CHECK( !container1.equal( container2 ));
+    }
+}
+

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 <http://www.gnu.org/licenses/>.
+ */
+
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE FLOAT_CMP_TESTS
+#include <boost/test/unit_test.hpp>
+
+#include <stdexcept>
+
+#include <opm/common/util/numeric/cmp.hpp>
+
+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 ));
+    }
+}
+
+