diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index 18dd617f..f37a71a9 100644 --- a/CMakeLists_files.cmake +++ b/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 diff --git a/opm/core/io/eclipse/EclipseIOUtil.hpp b/opm/core/io/eclipse/EclipseIOUtil.hpp new file mode 100644 index 00000000..1b8d3c20 --- /dev/null +++ b/opm/core/io/eclipse/EclipseIOUtil.hpp @@ -0,0 +1,34 @@ +#ifndef ECLIPSE_IO_UTIL_HPP +#define ECLIPSE_IO_UTIL_HPP + +#include +#include +#include + +namespace Opm +{ +namespace EclipseIOUtil +{ + + template + void addToStripedData(const std::vector& data, std::vector& 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 + void extractFromStripedData(const std::vector& data, std::vector& 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 diff --git a/opm/core/io/eclipse/EclipseReader.cpp b/opm/core/io/eclipse/EclipseReader.cpp index 9e122c37..74a2903d 100644 --- a/opm/core/io/eclipse/EclipseReader.cpp +++ b/opm/core/io/eclipse/EclipseReader.cpp @@ -19,6 +19,11 @@ #include "EclipseReader.hpp" #include +#include +#include +#include +#include + #include #include @@ -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 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 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 diff --git a/opm/core/io/eclipse/EclipseReader.hpp b/opm/core/io/eclipse/EclipseReader.hpp index 4b4932c9..7df97404 100644 --- a/opm/core/io/eclipse/EclipseReader.hpp +++ b/opm/core/io/eclipse/EclipseReader.hpp @@ -3,6 +3,9 @@ #include #include +#include +#include +#include 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 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 diff --git a/tests/test_readWriteWellStateData.cpp b/tests/test_readWriteWellStateData.cpp index c9dbc01d..791dd56c 100644 --- a/tests/test_readWriteWellStateData.cpp +++ b/tests/test_readWriteWellStateData.cpp @@ -27,8 +27,10 @@ #include #include +#include #include #include +#include #include #include #include @@ -251,6 +253,7 @@ BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData) simTimer->init(eclipseState->getSchedule()->getTimeMap()); eclipseWriter->writeInit(*simTimer); std::shared_ptr 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 blackoilState = createBlackOilState(eclipseState->getEclipseGrid()); wellState->init(wells, *blackoilState); + //Set test data for pressure + std::vector& pressure = blackoilState->pressure(); + for (std::vector::iterator iter = pressure.begin(); iter != pressure.end(); ++iter) { + *iter = 6.0; + } + + //Set test data for temperature + std::vector& temperature = blackoilState->temperature(); + for (std::vector::iterator iter = temperature.begin(); iter != temperature.end(); ++iter) { + *iter = 7.0; + } + + //Set test data for saturation water + std::vector 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 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 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 blackoilStateRestored = createBlackOilState(eclipseState->getEclipseGrid()); + Opm::restoreSOLUTIONData(restart_filename, 1, *eclipseState, *gridManager.c_grid(), phaseUsage, *blackoilStateRestored); + + std::vector swat_restored; + std::vector swat; + std::vector sgas_restored; + std::vector 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); }