Merge pull request #854 from chflo/OPM-215

Opm 215 EclipseReader restart file read of solution data

CMakeLists_files.cmake

@@ -286,6 +286,7 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/grid/cpgpreprocess/preprocess.h
 	opm/core/grid/cpgpreprocess/uniquepoints.h
 	opm/core/io/eclipse/CornerpointChopper.hpp
+	opm/core/io/eclipse/EclipseIOUtil.hpp
 	opm/core/io/eclipse/EclipseGridInspector.hpp
 	opm/core/io/eclipse/EclipseUnits.hpp
 	opm/core/io/eclipse/EclipseWriter.hpp

opm/core/io/eclipse/EclipseIOUtil.hpp

@@ -0,0 +1,34 @@
+#ifndef ECLIPSE_IO_UTIL_HPP
+#define ECLIPSE_IO_UTIL_HPP
+
+#include <vector>
+#include <string>
+#include <iostream>
+
+namespace Opm
+{
+namespace EclipseIOUtil
+{
+
+    template <typename T>
+    void addToStripedData(const std::vector<T>& data, std::vector<T>& result, size_t offset, size_t stride) {
+        int dataindex = 0;
+        for (size_t index = offset; index < result.size(); index += stride) {
+            result[index] = data[dataindex];
+            ++dataindex;
+        }
+    }
+
+
+    template <typename T>
+    void extractFromStripedData(const std::vector<T>& data, std::vector<T>& result, size_t offset, size_t stride) {
+        for (size_t index = offset; index < data.size(); index += stride) {
+            result.push_back(data[index]);
+        }
+    }
+
+
+} //namespace EclipseIOUtil
+} //namespace Opm
+
+#endif //ECLIPSE_IO_UTIL_HPP

opm/core/io/eclipse/EclipseReader.cpp

@@ -19,6 +19,11 @@
 
 #include "EclipseReader.hpp"
 #include <opm/core/simulator/WellState.hpp>
+#include <opm/core/simulator/SimulatorState.hpp>
+#include <opm/core/utility/Units.hpp>
+#include <opm/core/grid/GridHelpers.hpp>
+#include <opm/core/io/eclipse/EclipseIOUtil.hpp>
+
 #include <algorithm>
 
 #include <ert/ecl/ecl_file.h>
@@ -48,4 +53,105 @@ namespace Opm
             ecl_file_close(file_type);
         }
     }
+
+
+
+namespace {
+    void restoreTemperatureData(const ecl_file_type* file,
+                                const EclipseState& eclipse_state,
+                                const UnstructuredGrid& grid,
+                                SimulatorState& simulator_state) {
+        const char* temperature_keyword = "TEMP";
+
+        if (ecl_file_has_kw(file , temperature_keyword)) {
+            ecl_kw_type* temperature_kw = ecl_file_iget_named_kw(file, temperature_keyword, 0);
+
+            if (ecl_kw_get_size(temperature_kw) != Opm::UgGridHelpers::numCells(grid)) {
+                throw std::runtime_error("Read of restart file: Could not restore temperature data, length of data from file not equal number of cells");
+            }
+
+            float* temperature_data = ecl_kw_get_float_ptr(temperature_kw);
+
+            // factor and offset from the temperature values given in the deck to Kelvin
+            double scaling = eclipse_state.getDeckUnitSystem()->parse("Temperature")->getSIScaling();
+            double offset  = eclipse_state.getDeckUnitSystem()->parse("Temperature")->getSIOffset();
+
+            for (size_t index = 0; index < simulator_state.temperature().size(); ++index) {
+                simulator_state.temperature()[index] = unit::convert::from((double)temperature_data[index] - offset, scaling);
+            }
+        }
+    }
+
+
+    void restorePressureData(const ecl_file_type* file,
+                             const EclipseState& eclipse_state,
+                             const UnstructuredGrid& grid,
+                             SimulatorState& simulator_state) {
+        const char* pressure_keyword = "PRESSURE";
+
+        if (ecl_file_has_kw(file , pressure_keyword)) {
+
+            ecl_kw_type* pressure_kw = ecl_file_iget_named_kw(file, pressure_keyword, 0);
+
+            if (ecl_kw_get_size(pressure_kw) != Opm::UgGridHelpers::numCells(grid)) {
+                throw std::runtime_error("Read of restart file: Could not restore pressure data, length of data from file not equal number of cells");
+            }
+
+            float* pressure_data = ecl_kw_get_float_ptr(pressure_kw);
+            const double deck_pressure_unit = (eclipse_state.getDeckUnitSystem()->getType() == UnitSystem::UNIT_TYPE_METRIC) ? Opm::unit::barsa : Opm::unit::psia;
+            for (size_t index = 0; index < simulator_state.pressure().size(); ++index) {
+                simulator_state.pressure()[index] = unit::convert::from((double)pressure_data[index], deck_pressure_unit);
+            }
+        }
+    }
+}
+
+
+    void restoreSOLUTIONData(const std::string& restart_filename,
+                             int reportstep,
+                             const EclipseState& eclipseState,
+                             const UnstructuredGrid& grid,
+                             const PhaseUsage& phaseUsage,
+                             SimulatorState& simulator_state)
+    {
+        const char* filename = restart_filename.c_str();
+        ecl_file_type* file_type = ecl_file_open(filename, 0);
+
+        if (file_type != NULL) {
+            bool block_selected = ecl_file_select_rstblock_report_step(file_type , reportstep);
+
+            if (block_selected) {
+
+                restorePressureData(file_type, eclipseState, grid, simulator_state);
+                restoreTemperatureData(file_type, eclipseState, grid, simulator_state);
+
+                int numcells = Opm::UgGridHelpers::numCells(grid);
+
+                float* sgas_data = NULL;
+                float* swat_data = NULL;
+
+                if (phaseUsage.phase_used[BlackoilPhases::Aqua]) {
+                    ecl_kw_type* swat_kw = ecl_file_iget_named_kw(file_type , "SWAT", 0);
+                    swat_data = ecl_kw_get_float_ptr(swat_kw);
+                    std::vector<double> swat_datavec(&swat_data[0], &swat_data[numcells]);
+                    EclipseIOUtil::addToStripedData(swat_datavec, simulator_state.saturation(), phaseUsage.phase_pos[BlackoilPhases::Aqua], phaseUsage.num_phases);
+                }
+
+                if (phaseUsage.phase_used[BlackoilPhases::Vapour]) {
+                    ecl_kw_type* sgas_kw = ecl_file_iget_named_kw(file_type , "SGAS", 0);
+                    sgas_data = ecl_kw_get_float_ptr(sgas_kw);
+                    std::vector<double> sgas_datavec(&sgas_data[0], &sgas_data[numcells]);
+                    EclipseIOUtil::addToStripedData(sgas_datavec, simulator_state.saturation(), phaseUsage.phase_pos[BlackoilPhases::Vapour], phaseUsage.num_phases);
+                }
+            } else {
+                std::cerr << "Warning: Could not load solution data for report step " << reportstep << ", data for reportstep not found on file " << filename << std::endl;
+            }
+
+            ecl_file_close(file_type);
+        }
+    }
+
+
+
+
 } // namespace Opm

opm/core/io/eclipse/EclipseReader.hpp

@@ -3,6 +3,9 @@
 
 #include <string>
 #include <opm/core/simulator/WellState.hpp>
+#include <opm/core/simulator/SimulatorState.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
+#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
 
 namespace Opm
 {
@@ -18,6 +21,14 @@ namespace Opm
 ///     An instance of a WellState object, with correct size for each of the 5 contained std::vector<double> objects.
 ///
     void restoreOPM_XWELKeyword(const std::string& restart_filename, int report_step, WellState& wellState);
+    void restoreSOLUTIONData(const std::string& restart_filename,
+                             int report_step,
+                             const EclipseState &eclipseState,
+                             const UnstructuredGrid& grid,
+                             const PhaseUsage& phaseUsage,
+                             SimulatorState& simulator_state);
+
+
 }
 
 #endif // ECLIPSEREADER_HPP

tests/test_readWriteWellStateData.cpp

@@ -27,8 +27,10 @@
 
 #include <opm/core/io/eclipse/EclipseWriter.hpp>
 #include <opm/core/io/eclipse/EclipseReader.hpp>
+#include <opm/core/io/eclipse/EclipseIOUtil.hpp>
 #include <opm/core/grid/GridManager.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
+#include <opm/core/props/BlackoilPhases.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/simulator/SimulatorTimer.hpp>
@@ -251,6 +253,7 @@ BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData)
     simTimer->init(eclipseState->getSchedule()->getTimeMap());
     eclipseWriter->writeInit(*simTimer);
     std::shared_ptr<Opm::WellState> wellState(new Opm::WellState());
+    Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck);
 
     Opm::GridManager gridManager(deck);
     Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
@@ -258,6 +261,26 @@ BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData)
     std::shared_ptr<Opm::BlackoilState> blackoilState = createBlackOilState(eclipseState->getEclipseGrid());
     wellState->init(wells, *blackoilState);
 
+    //Set test data for pressure
+    std::vector<double>& pressure = blackoilState->pressure();
+    for (std::vector<double>::iterator iter = pressure.begin(); iter != pressure.end(); ++iter) {
+        *iter = 6.0;
+    }
+
+    //Set test data for temperature
+    std::vector<double>& temperature = blackoilState->temperature();
+    for (std::vector<double>::iterator iter = temperature.begin(); iter != temperature.end(); ++iter) {
+        *iter = 7.0;
+    }
+
+    //Set test data for saturation water
+    std::vector<double> swatdata(1000, 8);
+    Opm::EclipseIOUtil::addToStripedData(swatdata, blackoilState->saturation(), phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua], phaseUsage.num_phases);
+
+    //Set test data for saturation gas
+    std::vector<double> sgasdata(1000, 9);
+    Opm::EclipseIOUtil::addToStripedData(sgasdata, blackoilState->saturation(), phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour], phaseUsage.num_phases);
+
     setValuesInWellState(wellState);
     simTimer->setCurrentStepNum(1);
     eclipseWriter->writeTimeStep(*simTimer, *blackoilState, *wellState , false);
@@ -267,6 +290,8 @@ BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData)
     const std::string restart_filename = test_area_path + slash + eclipse_restart_filename;
     std::shared_ptr<Opm::WellState> wellStateRestored(new Opm::WellState());
     wellStateRestored->init(wells, *blackoilState);
+
+    //Read and verify OPM XWEL data
     Opm::restoreOPM_XWELKeyword(restart_filename, 1, *wellStateRestored);
 
     BOOST_CHECK_EQUAL_COLLECTIONS(wellState->bhp().begin(), wellState->bhp().end(), wellStateRestored->bhp().begin(), wellStateRestored->bhp().end());
@@ -275,5 +300,26 @@ BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData)
     BOOST_CHECK_EQUAL_COLLECTIONS(wellState->perfRates().begin(), wellState->perfRates().end(), wellStateRestored->perfRates().begin(), wellStateRestored->perfRates().end());
     BOOST_CHECK_EQUAL_COLLECTIONS(wellState->perfPress().begin(), wellState->perfPress().end(), wellStateRestored->perfPress().begin(), wellStateRestored->perfPress().end());
 
+
+    //Read and verify pressure, temperature and saturation data
+    std::shared_ptr<Opm::BlackoilState> blackoilStateRestored = createBlackOilState(eclipseState->getEclipseGrid());
+    Opm::restoreSOLUTIONData(restart_filename, 1, *eclipseState, *gridManager.c_grid(), phaseUsage, *blackoilStateRestored);
+
+    std::vector<double> swat_restored;
+    std::vector<double> swat;
+    std::vector<double> sgas_restored;
+    std::vector<double> sgas;
+    Opm::EclipseIOUtil::extractFromStripedData(blackoilStateRestored->saturation(), swat_restored, phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua], phaseUsage.num_phases);
+    Opm::EclipseIOUtil::extractFromStripedData(blackoilState->saturation(), swat, phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua], phaseUsage.num_phases);
+    Opm::EclipseIOUtil::extractFromStripedData(blackoilStateRestored->saturation(), sgas_restored, phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour], phaseUsage.num_phases);
+    Opm::EclipseIOUtil::extractFromStripedData(blackoilState->saturation(), sgas, phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour], phaseUsage.num_phases);
+
+    for (size_t cellindex = 0; cellindex < 1000; ++cellindex) {
+        BOOST_CHECK_CLOSE(blackoilState->pressure()[cellindex], blackoilStateRestored->pressure()[cellindex], 0.00001);
+        BOOST_CHECK_CLOSE(blackoilState->temperature()[cellindex], blackoilStateRestored->temperature()[cellindex], 0.00001);
+        BOOST_CHECK_CLOSE(swat[cellindex], swat_restored[cellindex], 0.00001);
+        BOOST_CHECK_CLOSE(sgas[cellindex], sgas_restored[cellindex], 0.00001);
+    }
+
     test_work_area_free(test_area);
 }