diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index f7e6c25c..e2ed1b1b 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -37,6 +37,7 @@ list (APPEND MAIN_SOURCE_FILES opm/core/grid/cpgpreprocess/geometry.c opm/core/grid/cpgpreprocess/preprocess.c opm/core/grid/cpgpreprocess/uniquepoints.c + opm/core/io/eclipse/BlackoilEclipseOutputWriter.cpp opm/core/io/eclipse/EclipseGridInspector.cpp opm/core/io/eclipse/EclipseGridParser.cpp opm/core/io/eclipse/writeECLData.cpp @@ -224,6 +225,7 @@ list (APPEND PUBLIC_HEADER_FILES opm/core/grid/cpgpreprocess/grdecl.h opm/core/grid/cpgpreprocess/preprocess.h opm/core/grid/cpgpreprocess/uniquepoints.h + opm/core/io/eclipse/BlackoilEclipseOutputWriter.hpp opm/core/io/eclipse/CornerpointChopper.hpp opm/core/io/eclipse/EclipseGridInspector.hpp opm/core/io/eclipse/EclipseGridParser.hpp diff --git a/opm/core/io/eclipse/BlackoilEclipseOutputWriter.cpp b/opm/core/io/eclipse/BlackoilEclipseOutputWriter.cpp new file mode 100644 index 00000000..1931c04b --- /dev/null +++ b/opm/core/io/eclipse/BlackoilEclipseOutputWriter.cpp @@ -0,0 +1,267 @@ + /* + Copyright 2013 Andreas Lauser + Copyright 2013 SINTEF ICT, Applied Mathematics. + + 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 "config.h" + +#include "BlackoilEclipseOutputWriter.hpp" + +#include + +#include +#include +#include +#include +#include + +#ifdef HAVE_ERT +#include +#include +#include +#include +#include +#include +#include +#include +#endif + +namespace Opm { +void BlackoilEclipseOutputWriter::writeInitFile(const SimulatorTimer &timer) +{ +#if HAVE_ERT + int phases = ECL_OIL_PHASE + ECL_WATER_PHASE; + bool endian_flip = true;//ECL_ENDIAN_FLIP; + bool fmt_file = false; + ecl_file_enum file_type = ECL_EGRID_FILE; + + ecl_grid_type* ecl_grid = newEclGrid_(); + char* gridFileName = ecl_util_alloc_filename(outputDir_.c_str(), baseName_.c_str(), file_type, fmt_file, timer.currentStepNum()); + fortio_type* fortio; + ecl_grid_fwrite_EGRID(ecl_grid, gridFileName); + free(gridFileName); + + char* initFileName = ecl_util_alloc_filename(outputDir_.c_str(), baseName_.c_str(), /*file_type=*/ECL_INIT_FILE, fmt_file, timer.currentStepNum()); + if (!ecl_util_fmt_file(initFileName, &fmt_file)) { + OPM_THROW(std::runtime_error, + "Could not determine formatted/unformatted status of file:" << initFileName << " non-standard name?" << std::endl); + } + fortio = fortio_open_writer(initFileName, fmt_file, endian_flip); + { + ecl_kw_type* poro_kw = newEclKeyword_(PORO_KW, ECL_FLOAT_TYPE); + time_t start_date; + + { + boost::posix_time::ptime start_date_(timer.currentDateTime()); + + tm td_tm = boost::posix_time::to_tm(start_date_); + start_date = mktime(&td_tm); + } + + ecl_init_file_fwrite_header(fortio, ecl_grid, poro_kw, phases, start_date); + ecl_kw_free(poro_kw); + } + + /* This collection of keywords is somewhat arbitrary and random. */ + saveEclKeyword_(fortio, "PERMX", ECL_FLOAT_TYPE); + saveEclKeyword_(fortio, "PERMY", ECL_FLOAT_TYPE); + saveEclKeyword_(fortio, "PERMZ", ECL_FLOAT_TYPE); + + fortio_fclose(fortio); + free(initFileName); +#else + OPM_THROW(std::runtime_error, + "The ERT libraries are required to write ECLIPSE output files."); +#endif // HAVE_ERT +} + +void BlackoilEclipseOutputWriter::writeReservoirState(const BlackoilState& reservoirState, const SimulatorTimer& timer) +{ +#if HAVE_ERT + ecl_file_enum file_type = ECL_UNIFIED_RESTART_FILE; // Alternatively ECL_RESTART_FILE for multiple restart files. + bool fmt_file = false; + + char *fileName = ecl_util_alloc_filename(outputDir_.c_str(), + baseName_.c_str(), + /*file_type=*/ECL_UNIFIED_RESTART_FILE, + fmt_file, + timer.currentStepNum()); + std::cout << "writing: " << fileName << "\n"; + int phases = ECL_OIL_PHASE + ECL_GAS_PHASE + ECL_WATER_PHASE; + double days = Opm::unit::convert::to(timer.currentTime(), Opm::unit::day); + time_t date = 0; + int nx = grid_.cartdims[0]; + int ny = grid_.cartdims[1]; + int nz = grid_.cartdims[2]; + int nactive = grid_.number_of_cells; + ecl_rst_file_type* rst_file; + + { + using namespace boost::posix_time; + ptime t0(boost::gregorian::date(1970, 1, 1)); + time_duration::sec_type seconds = (timer.currentDateTime() - t0).total_seconds(); + + date = time_t(seconds); + } + + if (timer.currentStepNum() > 0 && file_type == ECL_UNIFIED_RESTART_FILE) + rst_file = ecl_rst_file_open_append(fileName); + else + rst_file = ecl_rst_file_open_write(fileName); + + ecl_rst_file_fwrite_header(rst_file, timer.currentStepNum(), date, days, nx, ny, nz, nactive, phases); + ecl_rst_file_start_solution(rst_file); + + { + ecl_kw_type* pressure_kw = eclKeywordWrapper_("PRESSURE", reservoirState.pressure(), + /*offset=*/0, /*stride=*/1); + ecl_rst_file_add_kw(rst_file, pressure_kw); + ecl_kw_free(pressure_kw); + } + + { + ecl_kw_type* swat_kw = eclKeywordWrapper_("SWAT", reservoirState.saturation(), /*offset=*/0, /*stride=*/3); + ecl_rst_file_add_kw(rst_file, swat_kw); + ecl_kw_free(swat_kw); + } + + { + ecl_kw_type* soil_kw = eclKeywordWrapper_("SOIL", reservoirState.saturation(), /*offset=*/1, /*stride=*/3); + ecl_rst_file_add_kw(rst_file, soil_kw); + ecl_kw_free(soil_kw); + } + + { + ecl_kw_type* sgas_kw = eclKeywordWrapper_("SGAS", reservoirState.saturation(), /*offset=*/2, /*stride=*/3); + ecl_rst_file_add_kw(rst_file, sgas_kw); + ecl_kw_free(sgas_kw); + } + + ecl_rst_file_end_solution(rst_file); + ecl_rst_file_close(rst_file); + free(fileName); +#else + OPM_THROW(std::runtime_error, + "The ERT libraries are required to write ECLIPSE output files."); +#endif // HAVE_ERT +} + +void BlackoilEclipseOutputWriter::writeWellState(const WellState& wellState) +{ + for (int i = 0; i < wellState.perfPress().size(); ++i) + std::cout << "Perforation pressure " << i << ": " << wellState.perfPress()[i] << "\n"; + for (int i = 0; i < wellState.perfRates().size(); ++i) + std::cout << "Well rate " << i << ": " << wellState.perfRates()[i] << "\n"; +} + +#if HAVE_ERT +ecl_grid_type* BlackoilEclipseOutputWriter::newEclGrid_() { + if (eclipseParser_.hasField("DXV")) { + // make sure that the DYV and DZV keywords are present if the + // DXV keyword is used in the deck... + assert(eclipseParser_.hasField("DYV")); + assert(eclipseParser_.hasField("DZV")); + + const auto &dxv = eclipseParser_.getFloatingPointValue("DXV"); + const auto &dyv = eclipseParser_.getFloatingPointValue("DYV"); + const auto &dzv = eclipseParser_.getFloatingPointValue("DZV"); + + // creating a C array out of std::vector like this is pretty + // hacky and might even be unportable. having said that, it + // probably works with all currently known STL + // implementations... + return ecl_grid_alloc_dxv_dyv_dzv(dxv.size(), dyv.size(), dzv.size(), + &dxv[0], &dyv[0], &dzv[0], + /*actnum=*/NULL); + } + if (eclipseParser_.hasField("ZCORN")) { + struct grdecl grdecl = eclipseParser_.get_grdecl(); + + ecl_kw_type * coord_kw = newEclKeyword_(COORD_KW, ECL_FLOAT_TYPE); + ecl_kw_type * zcorn_kw = newEclKeyword_(ZCORN_KW, ECL_FLOAT_TYPE); + ecl_kw_type * actnum_kw = newEclKeyword_(ACTNUM_KW, ECL_INT_TYPE); + ecl_kw_type * mapaxes_kw = NULL; + + ecl_grid_type * grid ; + if (grdecl.mapaxes != NULL) + mapaxes_kw = newEclKeyword_(MAPAXES_KW, ECL_FLOAT_TYPE); + + grid = ecl_grid_alloc_GRDECL_kw(grdecl.dims[0], grdecl.dims[1], grdecl.dims[2], zcorn_kw, coord_kw, actnum_kw, mapaxes_kw); + + ecl_kw_free(coord_kw); + ecl_kw_free(zcorn_kw); + ecl_kw_free(actnum_kw); + if (mapaxes_kw != NULL) + ecl_kw_free(mapaxes_kw); + + return grid; + } + OPM_THROW(std::runtime_error, + "Can't create an ERT grid (no supported keywords found in deck)"); +} + +ecl_kw_type *BlackoilEclipseOutputWriter::newEclKeyword_(const std::string& keyword, ecl_type_enum ecl_type) const +{ + ecl_kw_type* ecl_kw = NULL; + + if (ecl_type == ECL_INT_TYPE) { + std::vector data = eclipseParser_.getIntegerValue(keyword); + ecl_kw = ecl_kw_alloc(keyword.c_str(), data.size(), ecl_type); + ecl_kw_set_memcpy_data(ecl_kw, &data[0]); + } else { + std::vector data = eclipseParser_.getFloatingPointValue(keyword); + if (ecl_type == ECL_DOUBLE_TYPE) { + ecl_kw = ecl_kw_alloc(keyword.c_str(), data.size(), ecl_type); + ecl_kw_set_memcpy_data(ecl_kw, &data[0]); + } else if (ecl_type == ECL_FLOAT_TYPE) { + ecl_kw = ecl_kw_alloc(keyword.c_str(), data.size(), ecl_type); + for (std::vector::size_type i=0; i < data.size(); i++) + ecl_kw_iset_float(ecl_kw, i, data[i]); + } + } + return ecl_kw; +} + +// TODO (?): unify with newEclKeyword_() +ecl_kw_type* BlackoilEclipseOutputWriter::eclKeywordWrapper_(const std::string& kw_name, + const std::vector &data, + int offset, + int stride) +{ + if (stride <= 0) + OPM_THROW(std::runtime_error, "Vector strides must be positive. Got stride = " << stride); + if ((stride * std::vector::size_type(grid_.number_of_cells)) != data.size()) + OPM_THROW(std::runtime_error, "Internal mismatch grid_.number_of_cells: " << grid_.number_of_cells << " data size: " << data.size() / stride); + + ecl_kw_type * ecl_kw = ecl_kw_alloc(kw_name.c_str(), grid_.number_of_cells, ECL_FLOAT_TYPE); + for (int i=0; i < grid_.number_of_cells; i++) + ecl_kw_iset_float(ecl_kw, i, data[i*stride + offset]); + return ecl_kw; +} + +void BlackoilEclipseOutputWriter::saveEclKeyword_(fortio_type* fortio, const std::string& kw, ecl_type_enum ecl_type) +{ + ecl_kw_type * ecl_kw = newEclKeyword_(kw, ecl_type); + if (ecl_kw != NULL) { + ecl_kw_fwrite(ecl_kw, fortio); + ecl_kw_free(ecl_kw); + } +} + +#endif // HAVE_ERT + +} // namespace Opm diff --git a/opm/core/io/eclipse/BlackoilEclipseOutputWriter.hpp b/opm/core/io/eclipse/BlackoilEclipseOutputWriter.hpp new file mode 100644 index 00000000..f940bfbd --- /dev/null +++ b/opm/core/io/eclipse/BlackoilEclipseOutputWriter.hpp @@ -0,0 +1,109 @@ +/* + Copyright 2013 Andreas Lauser + + 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 OPM_BLACKOIL_ECLIPSE_OUTPUT_WRITER_HPP +#define OPM_BLACKOIL_ECLIPSE_OUTPUT_WRITER_HPP + +#include +#include +#include +#include +#include + +#ifdef HAVE_ERT +#include +#include +#include +#include +#include +#include +#include +#endif + +#include + +namespace Opm { +/*! + * \brief A class to write the reservoir state and the well state of a + * blackoil simulation to disk using the Eclipse binary format. + * + * This class only writes files if the 'write_output' parameter is set + * to 1. It needs the ERT libraries to write to disk, so if the + * 'write_output' parameter is set but ERT is not available, all + * methods throw a std::runtime_error. + */ +class BlackoilEclipseOutputWriter +{ +public: + /*! + * \brief Sets the common attributes required to write eclipse + * binary files using ERT. + */ + BlackoilEclipseOutputWriter(const EclipseGridParser& eclipseParser, + const UnstructuredGrid& grid, + const std::string &outputDir, + const std::string &baseName) + : eclipseParser_(eclipseParser) + , grid_(grid) + , outputDir_(outputDir) + , baseName_(baseName) + {} + + /*! + * \brief Write the static eclipse data (grid, PVT curves, etc) to disk + */ + void writeInitFile(const SimulatorTimer &timer); + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection with + * visualization tools like ResInsight + * + * \param[in] reservoirState The thermodynamic state of the reservoir + */ + void writeReservoirState(const BlackoilState& reservoirState, + const SimulatorTimer& timer); + + /*! + * \brief Write a well state to disk for later inspection with + * visualization tools + * + * \param[in] wellState The production/injection data for all wells + */ + void writeWellState(const WellState& wellState); + +private: + const EclipseGridParser& eclipseParser_; + const UnstructuredGrid& grid_; + std::string outputDir_; + std::string baseName_; + +#if HAVE_ERT + ecl_grid_type* newEclGrid_(); + ecl_kw_type* eclKeywordWrapper_(const std::string& kw_name, + const std::vector &data, + int offset, + int stride); + ecl_kw_type* newEclKeyword_(const std::string &keyword , ecl_type_enum ecl_type) const; + void saveEclKeyword_(fortio_type* fortio, const std::string& keyword, ecl_type_enum ecl_type); +#endif +}; +} // namespace Opm + + +#endif // OPM_BLACKOIL_ECLIPSE_OUTPUT_WRITER_HPP