Merge pull request #1587 from tskille/eclio_huge_arrays

Adding support for larger EclFile arrays
2020-03-17 09:12:00 +01:00
parent 0e8c5faba9 e3ae38feb9
commit c0e328fd65
9 changed files with 185 additions and 58 deletions
--- a/opm/io/eclipse/EclFile.hpp
+++ b/opm/io/eclipse/EclFile.hpp
@@ -52,7 +52,7 @@ public:
      char_array.clear();
    }

-    using EclEntry = std::tuple<std::string, eclArrType, int>;
+    using EclEntry = std::tuple<std::string, eclArrType, long int>;
    std::vector<EclEntry> getList() const;

    template <typename T>
@@ -78,7 +78,7 @@ protected:

    std::vector<std::string> array_name;
    std::vector<eclArrType> array_type;
-    std::vector<int> array_size;
+    std::vector<long int> array_size;

    std::vector<unsigned long int> ifStreamPos;

--- a/opm/io/eclipse/EclOutput.hpp
+++ b/opm/io/eclipse/EclOutput.hpp
@@ -26,6 +26,7 @@

 #include <opm/io/eclipse/EclIOdata.hpp>
 #include <opm/io/eclipse/PaddedOutputString.hpp>
+#include <iostream>

 namespace Opm { namespace EclIO { namespace OutputStream {
    class Restart;
@@ -78,7 +79,7 @@ public:
    friend class OutputStream::SummarySpecification;

 private:
-    void writeBinaryHeader(const std::string& arrName, int size, eclArrType arrType);
+    void writeBinaryHeader(const std::string& arrName, long int size, eclArrType arrType);

    template <typename T>
    void writeBinaryArray(const std::vector<T>& data);
@@ -94,6 +95,7 @@ private:
    void writeFormattedCharArray(const std::vector<std::string>& data);
    void writeFormattedCharArray(const std::vector<PaddedOutputString<8>>& data);

+    void writeArrayType(const eclArrType arrType);
    std::string make_real_string(float value) const;
    std::string make_doub_string(double value) const;

--- a/opm/io/eclipse/EclUtil.hpp
+++ b/opm/io/eclipse/EclUtil.hpp
@@ -27,6 +27,7 @@
 namespace Opm { namespace EclIO {

    int flipEndianInt(int num);
+    long int flipEndianLongInt(long int num);
    float flipEndianFloat(float num);
    double flipEndianDouble(double num);

--- a/src/opm/io/eclipse/EclFile.cpp
+++ b/src/opm/io/eclipse/EclFile.cpp
@@ -30,7 +30,9 @@
 #include <sstream>
 #include <string>
 #include <numeric>
+#include <cmath>

+#include <iostream>

 // anonymous namespace for EclFile

@@ -67,35 +69,60 @@ bool isEOF(std::fstream* fileH)
 }


-void readBinaryHeader(std::fstream& fileH, std::string& arrName,
-                      int& size, Opm::EclIO::eclArrType &arrType)
+void readBinaryHeader(std::fstream& fileH, std::string& tmpStrName,
+                      int& tmpSize, std::string& tmpStrType)
 {
    int bhead;
-    std::string tmpStrName(8,' ');
-    std::string tmpStrType(4,' ');

    fileH.read(reinterpret_cast<char*>(&bhead), sizeof(bhead));
    bhead = Opm::EclIO::flipEndianInt(bhead);

-    if (bhead != 16) {
+    if (bhead != 16){
        std::string message="Error reading binary header. Expected 16 bytes of header data, found " + std::to_string(bhead);
        OPM_THROW(std::runtime_error, message);
    }

    fileH.read(&tmpStrName[0], 8);

-    fileH.read(reinterpret_cast<char*>(&size), sizeof(size));
-    size = Opm::EclIO::flipEndianInt(size);
+    fileH.read(reinterpret_cast<char*>(&tmpSize), sizeof(tmpSize));
+    tmpSize = Opm::EclIO::flipEndianInt(tmpSize);

    fileH.read(&tmpStrType[0], 4);

    fileH.read(reinterpret_cast<char*>(&bhead), sizeof(bhead));
    bhead = Opm::EclIO::flipEndianInt(bhead);

-    if (bhead != 16) {
+    if (bhead != 16){
        std::string message="Error reading binary header. Expected 16 bytes of header data, found " + std::to_string(bhead);
        OPM_THROW(std::runtime_error, message);
    }
+}
+
+void readBinaryHeader(std::fstream& fileH, std::string& arrName,
+                      long int& size, Opm::EclIO::eclArrType &arrType)
+{
+    std::string tmpStrName(8,' ');
+    std::string tmpStrType(4,' ');
+    int tmpSize;
+
+    readBinaryHeader(fileH, tmpStrName, tmpSize, tmpStrType);
+    
+    if (tmpStrType == "X231"){
+        std::string x231ArrayName = tmpStrName;
+        int x231exp = tmpSize * (-1);
+
+        readBinaryHeader(fileH, tmpStrName, tmpSize, tmpStrType);
+        
+        if (x231ArrayName != tmpStrName)
+            OPM_THROW(std::runtime_error, "Invalied X231 header, name should be same in both headers'");
+
+        if (x231exp < 0)
+            OPM_THROW(std::runtime_error, "Invalied X231 header, size of array should be negative'");
+        
+        size = static_cast<long int>(tmpSize) + static_cast<long int>(x231exp) * pow(2,31); 
+    } else {
+        size = static_cast<long int>(tmpSize);
+    }

    arrName = tmpStrName;
    if (tmpStrType == "INTE")
@@ -114,8 +141,7 @@ void readBinaryHeader(std::fstream& fileH, std::string& arrName,
        OPM_THROW(std::runtime_error, "Error, unknown array type '" + tmpStrType +"'");
 }

-
-unsigned long int sizeOnDiskBinary(int num, Opm::EclIO::eclArrType arrType)
+unsigned long int sizeOnDiskBinary(long int num, Opm::EclIO::eclArrType arrType)
 {
    unsigned long int size = 0;

@@ -125,24 +151,32 @@ unsigned long int sizeOnDiskBinary(int num, Opm::EclIO::eclArrType arrType)
            OPM_THROW(std::invalid_argument, message);
        }
    } else {
-        auto sizeData = Opm::EclIO::block_size_data_binary(arrType);
-
-        int sizeOfElement = std::get<0>(sizeData);
-        int maxBlockSize = std::get<1>(sizeData);
-        int maxNumberOfElements = maxBlockSize / sizeOfElement;
-
-        size = num * sizeOfElement;
-        size = size + ((num-1) / maxNumberOfElements) * 2 * Opm::EclIO::sizeOfInte; // 8 byte (two integers) every 1000 element
-
        if (num > 0) {
-            size = size + 2 * Opm::EclIO::sizeOfInte;
+            auto sizeData = Opm::EclIO::block_size_data_binary(arrType);
+
+            int sizeOfElement = std::get<0>(sizeData);
+            int maxBlockSize = std::get<1>(sizeData);
+            int maxNumberOfElements = maxBlockSize / sizeOfElement;
+
+            long unsigned int numBlocks = static_cast<unsigned long int>(num)/static_cast<unsigned long int>(maxNumberOfElements);
+            long unsigned int rest = static_cast<unsigned long int>(num) - numBlocks*static_cast<unsigned long int>(maxNumberOfElements);
+
+            long unsigned int size2Inte = static_cast<long unsigned int>(Opm::EclIO::sizeOfInte) * 2;
+            long unsigned int sizeFullBlocks = numBlocks * (static_cast<long unsigned int>(maxBlockSize) + size2Inte);
+
+            long unsigned int sizeLastBlock = 0;
+
+            if (rest > 0)
+                sizeLastBlock = rest * static_cast<long unsigned int>(sizeOfElement) + size2Inte;
+
+            size = sizeFullBlocks + sizeLastBlock;
        }
    }

    return size;
 }

-unsigned long int sizeOnDiskFormatted(const int num, Opm::EclIO::eclArrType arrType)
+unsigned long int sizeOnDiskFormatted(const long int num, Opm::EclIO::eclArrType arrType)
 {
    unsigned long int size = 0;

@@ -189,7 +223,7 @@ unsigned long int sizeOnDiskFormatted(const int num, Opm::EclIO::eclArrType arrT


 template<typename T, typename T2>
-std::vector<T> readBinaryArray(std::fstream& fileH, const int size, Opm::EclIO::eclArrType type,
+std::vector<T> readBinaryArray(std::fstream& fileH, const long int size, Opm::EclIO::eclArrType type,
                               std::function<T(T2)>& flip)
 {
    std::vector<T> arr;
@@ -201,7 +235,7 @@ std::vector<T> readBinaryArray(std::fstream& fileH, const int size, Opm::EclIO::

    arr.reserve(size);

-    int rest = size;
+    long int rest = size;
    while (rest > 0) {
        int dhead;
        fileH.read(reinterpret_cast<char*>(&dhead), sizeof(dhead));
@@ -240,27 +274,27 @@ std::vector<T> readBinaryArray(std::fstream& fileH, const int size, Opm::EclIO::
 }


-std::vector<int> readBinaryInteArray(std::fstream &fileH, const int size)
+std::vector<int> readBinaryInteArray(std::fstream &fileH, const long int size)
 {
    std::function<int(int)> f = Opm::EclIO::flipEndianInt;
    return readBinaryArray<int,int>(fileH, size, Opm::EclIO::INTE, f);
 }


-std::vector<float> readBinaryRealArray(std::fstream& fileH, const int size)
+std::vector<float> readBinaryRealArray(std::fstream& fileH, const long int size)
 {
    std::function<float(float)> f = Opm::EclIO::flipEndianFloat;
    return readBinaryArray<float,float>(fileH, size, Opm::EclIO::REAL, f);
 }


-std::vector<double> readBinaryDoubArray(std::fstream& fileH, const int size)
+std::vector<double> readBinaryDoubArray(std::fstream& fileH, const long int size)
 {
    std::function<double(double)> f = Opm::EclIO::flipEndianDouble;
    return readBinaryArray<double,double>(fileH, size, Opm::EclIO::DOUB, f);
 }

-std::vector<bool> readBinaryLogiArray(std::fstream &fileH, const int size)
+std::vector<bool> readBinaryLogiArray(std::fstream &fileH, const long int size)
 {
    std::function<bool(unsigned int)> f = [](unsigned int intVal)
                                          {
@@ -279,7 +313,7 @@ std::vector<bool> readBinaryLogiArray(std::fstream &fileH, const int size)
 }


-std::vector<std::string> readBinaryCharArray(std::fstream& fileH, const int size)
+std::vector<std::string> readBinaryCharArray(std::fstream& fileH, const long int size)
 {
    using Char8 = std::array<char, 8>;
    std::function<std::string(Char8)> f = [](const Char8& val)
@@ -292,7 +326,7 @@ std::vector<std::string> readBinaryCharArray(std::fstream& fileH, const int size


 void readFormattedHeader(std::fstream& fileH, std::string& arrName,
-                         int &num, Opm::EclIO::eclArrType &arrType)
+                         long int &num, Opm::EclIO::eclArrType &arrType)
 {
    std::string line;
    std::getline(fileH,line);
@@ -310,7 +344,7 @@ void readFormattedHeader(std::fstream& fileH, std::string& arrName,
    std::string antStr = line.substr(p2 + 1, p3 - p2 - 1);
    std::string arrTypeStr = line.substr(p3 + 1, p4 - p3 - 1);

-    num = std::stoi(antStr);
+    num = std::stol(antStr);

    if (arrTypeStr == "INTE")
        arrType = Opm::EclIO::INTE;
@@ -357,7 +391,7 @@ std::vector<T> readFormattedArray(const std::string& file_str, const int size, l
 }


-std::vector<int> readFormattedInteArray(const std::string& file_str, const int size, long int fromPos)
+std::vector<int> readFormattedInteArray(const std::string& file_str, const long int size, long int fromPos)
 {

    std::function<int(const std::string&)> f = [](const std::string& val)
@@ -369,7 +403,7 @@ std::vector<int> readFormattedInteArray(const std::string& file_str, const int s
 }


-std::vector<std::string> readFormattedCharArray(const std::string& file_str, const int size, long int fromPos)
+std::vector<std::string> readFormattedCharArray(const std::string& file_str, const long int size, long int fromPos)
 {
    std::vector<std::string> arr;
    arr.reserve(size);
@@ -393,7 +427,7 @@ std::vector<std::string> readFormattedCharArray(const std::string& file_str, con
 }


-std::vector<float> readFormattedRealArray(const std::string& file_str, const int size, long int fromPos)
+std::vector<float> readFormattedRealArray(const std::string& file_str, const long int size, long int fromPos)
 {

    std::function<float(const std::string&)> f = [](const std::string& val)
@@ -408,7 +442,7 @@ std::vector<float> readFormattedRealArray(const std::string& file_str, const int
 }


-std::vector<bool> readFormattedLogiArray(const std::string& file_str, const int size, long int fromPos)
+std::vector<bool> readFormattedLogiArray(const std::string& file_str, const long int size, long int fromPos)
 {

    std::function<bool(const std::string&)> f = [](const std::string& val)
@@ -426,7 +460,7 @@ std::vector<bool> readFormattedLogiArray(const std::string& file_str, const int
    return readFormattedArray<bool>(file_str, size, fromPos, f);
 }

-std::vector<double> readFormattedDoubArray(const std::string& file_str, const int size, long int fromPos)
+std::vector<double> readFormattedDoubArray(const std::string& file_str, const long int size, long int fromPos)
 {

    std::function<double(const std::string&)> f = [](std::string val)
@@ -480,8 +514,8 @@ EclFile::EclFile(const std::string& filename, bool preload) : inputFilename(file
    while (!isEOF(&fileH)) {
        std::string arrName(8,' ');
        eclArrType arrType;
-        int num;
-
+        long int num;
+        
        if (formatted) {
            readFormattedHeader(fileH,arrName,num,arrType);
        } else {
@@ -499,14 +533,17 @@ EclFile::EclFile(const std::string& filename, bool preload) : inputFilename(file

        arrayLoaded.push_back(false);

-        if (formatted) {
-            unsigned long int sizeOfNextArray = sizeOnDiskFormatted(num, arrType);
-            fileH.seekg(static_cast<std::streamoff>(sizeOfNextArray), std::ios_base::cur);
-        } else {
-            unsigned long int sizeOfNextArray = sizeOnDiskBinary(num, arrType);
-            fileH.seekg(static_cast<std::streamoff>(sizeOfNextArray), std::ios_base::cur);
+        if (num > 0){ 
+            if (formatted) {
+                unsigned long int sizeOfNextArray = sizeOnDiskFormatted(num, arrType);
+                fileH.seekg(static_cast<std::streamoff>(sizeOfNextArray), std::ios_base::cur);
+            } else {
+                unsigned long int sizeOfNextArray = sizeOnDiskBinary(num, arrType);
+                fileH.seekg(static_cast<std::streamoff>(sizeOfNextArray), std::ios_base::cur);
+            }
        }

+
        n++;
    };

@@ -812,17 +849,26 @@ EclFile::seekPosition(const std::vector<std::string>::size_type arrIndex) const
    //
    //
    //   (*) unformatted header size = 24 bytes =
-    //            4 byte control character
+    //            4 byte => size of header in byte
    //          + 8 byte header/vector name
    //          + 4 byte for number of elements
    //          + 4 byte for element type
-    //          + 4 byte control character
+    //          + 4 byte size of header in byte
+    //
+    //      16 byte header (8+4+4)
    //
    //       +------+------------+------+------+------+
    //       | Ctrl | Keyword    | #elm | type | Ctrl |  (item)
    //       |  4   |  8         |  4   |  4   |  4   |  (#bytes)
    //       +------+------------+------+------+------+
    //
+    //      20 byte header (8+8+4)
+    //
+    //       +------+------------+------+------+------+
+    //       | Ctrl | Keyword    | #elm | type | Ctrl |  (item)
+    //       |  4   |  8         |  8   |  4   |  4   |  (#bytes)
+    //       +------+------------+------+------+------+
+    //

    const auto headerSize = this->formatted ? 30ul : 24ul;
    const auto datapos    = this->ifStreamPos[arrIndex];
--- a/src/opm/io/eclipse/EclOutput.cpp
+++ b/src/opm/io/eclipse/EclOutput.cpp
@@ -91,12 +91,28 @@ void EclOutput::flushStream()
    this->ofileH.flush();
 }

-void EclOutput::writeBinaryHeader(const std::string&arrName, int size, eclArrType arrType)
+void EclOutput::writeBinaryHeader(const std::string&arrName, long int size, eclArrType arrType)
 {
+    int bhead = flipEndianInt(16);
    std::string name = arrName + std::string(8 - arrName.size(),' ');

+    // write X231 header if size larger that limits for 4 byte integers
+    if (size > std::numeric_limits<int>::max()) {
+        long int val231 = std::pow(2,31);
+        long int x231 = size / val231;
+
+        int flippedx231 = flipEndianInt(static_cast<int>( (-1)*x231 ));
+        
+        ofileH.write(reinterpret_cast<char*>(&bhead), sizeof(bhead));
+        ofileH.write(name.c_str(), 8);
+        ofileH.write(reinterpret_cast<char*>(&flippedx231), sizeof(flippedx231));
+        ofileH.write("X231", 4);
+        ofileH.write(reinterpret_cast<char*>(&bhead), sizeof(bhead));
+        
+        size = size - (x231 * val231);
+    }
+
    int flippedSize = flipEndianInt(size);
-    int bhead = flipEndianInt(16);

    ofileH.write(reinterpret_cast<char*>(&bhead), sizeof(bhead));

@@ -127,18 +143,18 @@ void EclOutput::writeBinaryHeader(const std::string&arrName, int size, eclArrTyp
    ofileH.write(reinterpret_cast<char *>(&bhead), sizeof(bhead));
 }

-
 template <typename T>
 void EclOutput::writeBinaryArray(const std::vector<T>& data)
 {
-    int rest,num,rval;
+    int num, rval;
+    long int rest;
    int dhead;
    float value_f;
    double value_d;
    int intVal;

-    int n = 0;
-    int size = data.size();
+    long int n = 0;
+    long int size = data.size();

    eclArrType arrType = MESS;

@@ -162,13 +178,13 @@ void EclOutput::writeBinaryArray(const std::vector<T>& data)
        OPM_THROW(std::runtime_error, "fstream fileH not open for writing");
    }

-    rest = size * sizeOfElement;
+    rest = size * static_cast<long int>(sizeOfElement);
    while (rest > 0) {
        if (rest > maxBlockSize) {
            rest -= maxBlockSize;
            num = maxNumberOfElements;
        } else {
-            num = rest / sizeOfElement;
+            num = static_cast<int>(rest) / sizeOfElement;
            rest = 0;
        }

--- a/src/opm/io/eclipse/EclUtil.cpp
+++ b/src/opm/io/eclipse/EclUtil.cpp
@@ -30,6 +30,11 @@ int Opm::EclIO::flipEndianInt(int num)
    return static_cast<int>(tmp);
 }

+long int Opm::EclIO::flipEndianLongInt(long int num)
+{
+    unsigned long int tmp = __builtin_bswap64(num);
+    return static_cast<long int>(tmp);
+}

 float Opm::EclIO::flipEndianFloat(float num)
 {
--- a/test_util/convertECL.cpp
+++ b/test_util/convertECL.cpp
@@ -9,7 +9,7 @@
 #include <opm/io/eclipse/EclOutput.hpp>

 using namespace Opm::EclIO;
-using EclEntry = std::tuple<std::string, eclArrType, int>;
+using EclEntry = std::tuple<std::string, eclArrType, long int>;

 template <typename T>
 void write(EclOutput& outFile, EclFile& file1,
--- a/tests/test_ERst.cpp
+++ b/tests/test_ERst.cpp
@@ -123,7 +123,7 @@ BOOST_AUTO_TEST_CASE(TestERst_1) {


    // try to get a list of vectors from non-existing report step, should throw exception
-    std::vector<std::tuple<std::string, eclArrType, int>> rstArrays; // = rst1.listOfRstArrays(4);
+    std::vector<std::tuple<std::string, eclArrType, long int>> rstArrays; // = rst1.listOfRstArrays(4);
    BOOST_CHECK_THROW(rstArrays = rst1.listOfRstArrays(4), std::invalid_argument);

    // non exising report step number, should throw exception
--- a/tests/test_EclIO.cpp
+++ b/tests/test_EclIO.cpp
@@ -27,8 +27,10 @@
 #include <limits>
 #include <tuple>
 #include <cmath>
+#include <numeric>

 #include <opm/io/eclipse/EclFile.hpp>
+#include <opm/io/eclipse/EclUtil.hpp>
 #include "WorkArea.cpp"

 #include <opm/io/eclipse/EclOutput.hpp>
@@ -36,6 +38,8 @@
 #define BOOST_TEST_MODULE Test EclIO
 #include <boost/test/unit_test.hpp>

+
+
 using namespace Opm::EclIO;

 template<typename InputIterator1, typename InputIterator2>
@@ -71,6 +75,59 @@ bool operator==(const std::vector<T> & t1, const std::vector<T> & t2)
    return std::equal(t1.begin(), t1.end(), t2.begin(), t2.end());
 }

+void write_header(std::ofstream& ofileH, std::string& arrName, int size, std::string arrtype){
+
+    int bhead = flipEndianInt(16);
+    int fsize = flipEndianInt(size);
+
+    ofileH.write(reinterpret_cast<char*>(&bhead), sizeof(bhead));
+    ofileH.write(arrName.c_str(), 8);
+    ofileH.write(reinterpret_cast<char*>(&fsize), sizeof(fsize));
+    ofileH.write(arrtype.c_str(), 4);
+    ofileH.write(reinterpret_cast<char*>(&bhead), sizeof(bhead));
+}
+
+
+BOOST_AUTO_TEST_CASE(TestEclFile_X231) {
+
+    std::string filename = "TEST.DAT";
+    std::string arrName = "TESTX231";
+
+    std::vector<int> ivect(10);
+    std::iota(ivect.begin(), ivect.end(), -4);    
+
+    {
+        std::ofstream ofileH;
+        ofileH.open(filename, std::ios_base::binary);    
+
+        int size = static_cast<int>((-1) * std::pow(2,31) + 10); 
+
+        write_header(ofileH, arrName, -1, std::string("X231"));
+        write_header(ofileH, arrName, size, std::string("INTE"));
+
+        int sizeData = ivect.size()*sizeof(int);
+        sizeData = flipEndianInt(sizeData);
+        
+        ofileH.write(reinterpret_cast<char*>(&sizeData), sizeof(sizeData));
+
+        for (auto v : ivect){
+            int fval = flipEndianInt(v); 
+            ofileH.write(reinterpret_cast<char*>(&fval), sizeof(fval));
+        }
+        
+        ofileH.write(reinterpret_cast<char*>(&sizeData), sizeof(sizeData));
+        ofileH.close();
+    }
+    
+    EclFile test1(filename);
+    auto array = test1.get<int>(arrName);
+    
+    for (size_t n = 0; n < 10; n++){
+        BOOST_CHECK_EQUAL(array[n], ivect[n]);
+    }
+}
+
+
 BOOST_AUTO_TEST_CASE(TestEclFile_BINARY) {

    std::string testFile="ECLFILE.INIT";