diff --git a/opm/output/eclipse/EclipseWriter.hpp b/opm/output/eclipse/EclipseIO.hpp similarity index 67% rename from opm/output/eclipse/EclipseWriter.hpp rename to opm/output/eclipse/EclipseIO.hpp index 55b453539..409b0cac3 100644 --- a/opm/output/eclipse/EclipseWriter.hpp +++ b/opm/output/eclipse/EclipseIO.hpp @@ -42,13 +42,13 @@ class EclipseState; * \brief A class to write the reservoir state and the well state of a * blackoil simulation to disk using the Eclipse binary format. */ -class EclipseWriter { +class EclipseIO { public: /*! * \brief Sets the common attributes required to write eclipse * binary files using ERT. */ - EclipseWriter( const EclipseState&, EclipseGrid ); + EclipseIO( const EclipseState&, EclipseGrid ); @@ -84,23 +84,49 @@ public: * permeabilities KRO, KRW and KRG and fluxes. The keywords which * can be added here are represented with mnenonics in the RPTRST * keyword. - * - * By default the various solution fields like PRESSURE and - * SWAT/SGAS should be written in single precision (i.e. as - * float). That is what eclipse does, and probably what most third - * party application expect - however passing false for the - * optional variable write_float the solution fields will be - * written in double precision. */ + void writeTimeStep( int report_step, bool isSubstep, double seconds_elapsed, data::Solution, - data::Wells, - bool write_float = true); + data::Wells); - EclipseWriter( const EclipseWriter& ) = delete; - ~EclipseWriter(); + + /* + Will load solution data and wellstate from the restart + file. This method will consult the IOConfig object to get + filename and report step to restart from. + + The map keys should be a map of keyword names and their + corresponding dimension object, i.e. to load the state from a + simple two phase simulation you would pass: + + keys = {{"PRESSURE" , UnitSystem::measure::pressure}, + {"SWAT" , UnitSystem::measure::identity }} + + For a three phase black oil simulation you would add pairs for + SGAS, RS and RV. If you ask for keys which are not found in the + restart file an exception will be raised, the happens if the + size of a vector is wrong. + + The function will consult the InitConfig object in the + EclipseState object to determine which file and report step to + load. + + The return value is of type 'data::Solution', which is the same + container type which is used by the EclipseIO, but observe + that the dim and target elements carry no information: + + - The returned double data has been converted to SI. + . The target is unconditionally set to 'RESTART_SOLUTION' + */ + std::pair< data::Solution, data::Wells > + loadRestart(const std::map& keys) const; + + + EclipseIO( const EclipseIO& ) = delete; + ~EclipseIO(); private: class Impl; diff --git a/opm/output/eclipse/EclipseIOUtil.hpp b/opm/output/eclipse/EclipseIOUtil.hpp index 1b8d3c202..8f07d159b 100644 --- a/opm/output/eclipse/EclipseIOUtil.hpp +++ b/opm/output/eclipse/EclipseIOUtil.hpp @@ -5,6 +5,7 @@ #include #include + namespace Opm { namespace EclipseIOUtil diff --git a/opm/output/eclipse/EclipseReader.hpp b/opm/output/eclipse/EclipseReader.hpp deleted file mode 100644 index 8940399e9..000000000 --- a/opm/output/eclipse/EclipseReader.hpp +++ /dev/null @@ -1,34 +0,0 @@ -#ifndef ECLIPSEREADER_HPP -#define ECLIPSEREADER_HPP - -#include - -#include -#include -#include - -namespace Opm { - - class EclipseState; - -/// -/// \brief init_from_restart_file -/// Reading from the restart file, information stored under the OPM_XWEL keyword and SOLUTION data is in this method filled into -/// an instance of a wellstate object and a SimulatorState object. -/// \param grid -/// UnstructuredGrid reference -/// \param pu -/// PhaseUsage reference -/// \param simulator_state -/// An instance of a SimulatorState object -/// \param wellstate -/// An instance of a WellState object, with correct size for each of the 5 contained std::vector objects -/// - - std::pair< data::Solution, data::Wells > - init_from_restart_file( const EclipseState&, int numcells ); - - -} - -#endif // ECLIPSEREADER_HPP diff --git a/opm/output/eclipse/RestartIO.hpp b/opm/output/eclipse/RestartIO.hpp new file mode 100644 index 000000000..a2b19de2f --- /dev/null +++ b/opm/output/eclipse/RestartIO.hpp @@ -0,0 +1,91 @@ +/* + Copyright (c) 2016 Statoil ASA + Copyright (c) 2013-2015 Andreas Lauser + Copyright (c) 2013 SINTEF ICT, Applied Mathematics. + Copyright (c) 2013 Uni Research AS + Copyright (c) 2015 IRIS AS + + 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 RESTART_IO_HPP +#define RESTART_IO_HPP + +#include + +#include +#include +#include + +#include +#include +#include + +#include +#include +#include +#include + +namespace Opm { + +class EclipseGrid; +class EclipseState; +class Phases; +class Schedule; + +namespace RestartIO { + + +/* + The two loose functions RestartIO::save() and RestartIO::load() can + be used to save and load reservoir and well state from restart + files. Observe that these functions 'just do it', i.e. the checking + of which report step to load from, if output is enabled at all and + so on is handled by an outer scope. + + If the filename corresponds to unified eclipse restart file, + i.e. UNRST the functions will seek correctly to the correct report + step, and truncate in the case of save. For any other filename the + functions will start reading and writing from file offset zero. If + the input filename does not correspond to a unified restart file + there is no consistency checking between filename and report step; + i.e. these calls: + + load("CASE.X0010" , 99 , ...) + save("CASE.X0010" , 99 , ...) + + will read and write to the file "CASE.X0010" - completely ignoring + the report step argument '99'. +*/ + +void save(const std::string& filename, + int report_step, + double seconds_elapsed, + data::Solution cells, + data::Wells wells, + const EclipseState& es, + const EclipseGrid& grid); + + + +std::pair< data::Solution, data::Wells > load( const std::string& filename, + int report_step, + const std::map& keys, + const EclipseState& es, + const EclipseGrid& grid); + +} +} +#endif diff --git a/opm/output/eclipse/Summary.hpp b/opm/output/eclipse/Summary.hpp index f17677286..ad044cb42 100644 --- a/opm/output/eclipse/Summary.hpp +++ b/opm/output/eclipse/Summary.hpp @@ -43,8 +43,6 @@ namespace Opm { namespace out { - class RegionCache; - class Summary { public: Summary( const EclipseState&, const SummaryConfig&, const EclipseGrid& ); @@ -54,7 +52,6 @@ class Summary { void add_timestep( int report_step, double secs_elapsed, const EclipseState& es, - const RegionCache& regionCache, const data::Wells&, const data::Solution& ); @@ -67,6 +64,7 @@ class Summary { class keyword_handlers; const EclipseGrid& grid; + out::RegionCache regionCache; ERT::ert_unique_ptr< ecl_sum_type, ecl_sum_free > ecl_sum; std::unique_ptr< keyword_handlers > handlers; const ecl_sum_tstep_type* prev_tstep = nullptr; diff --git a/src/opm/output/eclipse/EclipseIO.cpp b/src/opm/output/eclipse/EclipseIO.cpp new file mode 100644 index 000000000..57a1d48ff --- /dev/null +++ b/src/opm/output/eclipse/EclipseIO.cpp @@ -0,0 +1,476 @@ +/* + Copyright (c) 2016 Statoil ASA + Copyright (c) 2013-2015 Andreas Lauser + Copyright (c) 2013 SINTEF ICT, Applied Mathematics. + Copyright (c) 2013 Uni Research AS + Copyright (c) 2015 IRIS AS + + 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 + +#include "config.h" + +#include "EclipseIO.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include // unique_ptr +#include // move + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define OPM_XWEL "OPM_XWEL" +#define OPM_IWEL "OPM_IWEL" + +// namespace start here since we don't want the ERT headers in it +namespace Opm { +namespace { + + + + +void writeKeyword( ERT::FortIO& fortio , + const std::string& keywordName, + const std::vector &data ) { + ERT::EclKW< int > kw( keywordName, data ); + kw.fwrite( fortio ); +} + +/* + This overload hardcodes the common assumption that properties which + are stored internally as double values in OPM should be stored as + float values in the ECLIPSE formatted binary files. +*/ + +void writeKeyword( ERT::FortIO& fortio , + const std::string& keywordName, + const std::vector &data) { + + ERT::EclKW< float > kw( keywordName, data ); + kw.fwrite( fortio ); + +} + + + + + +/// Convert OPM phase usage to ERT bitmask +inline int ertPhaseMask( const Phases& phase ) { + return ( phase.active( Phase::WATER ) ? ECL_WATER_PHASE : 0 ) + | ( phase.active( Phase::OIL ) ? ECL_OIL_PHASE : 0 ) + | ( phase.active( Phase::GAS ) ? ECL_GAS_PHASE : 0 ); +} + +class RFT { + public: + RFT( const std::string& output_dir, + const std::string& basename, + bool format ); + + void writeTimeStep( std::vector< const Well* >, + const EclipseGrid& grid, + int report_step, + time_t current_time, + double days, + const UnitSystem& units, + data::Solution cells); + private: + std::string filename; +}; + + + RFT::RFT( const std::string& output_dir, + const std::string& basename, + bool format ) : + filename( ERT::EclFilename( output_dir, basename, ECL_RFT_FILE, format ) ) +{} + + +void RFT::writeTimeStep( std::vector< const Well* > wells, + const EclipseGrid& grid, + int report_step, + time_t current_time, + double days, + const UnitSystem& units, + data::Solution cells) { + + using rft = ERT::ert_unique_ptr< ecl_rft_node_type, ecl_rft_node_free >; + const std::vector& pressure = cells.data("PRESSURE"); + const std::vector& swat = cells.data("SWAT"); + const std::vector& sgas = cells.has("SGAS") ? cells.data("SGAS") : std::vector( pressure.size() , 0 ); + ERT::FortIO fortio(filename, std::ios_base::out); + + cells.convertFromSI( units ); + for ( const auto& well : wells ) { + if( !( well->getRFTActive( report_step ) + || well->getPLTActive( report_step ) ) ) + continue; + + auto* rft_node = ecl_rft_node_alloc_new( well->name().c_str(), "RFT", + current_time, days ); + + for( const auto& completion : well->getCompletions( report_step ) ) { + const size_t i = size_t( completion.getI() ); + const size_t j = size_t( completion.getJ() ); + const size_t k = size_t( completion.getK() ); + + if( !grid.cellActive( i, j, k ) ) continue; + + const auto index = grid.activeIndex( i, j, k ); + const double depth = grid.getCellDepth( i, j, k ); + const double press = pressure[ index ]; + const double satwat = swat[ index ]; + const double satgas = sgas[ index ]; + + auto* cell = ecl_rft_cell_alloc_RFT( + i, j, k, depth, press, satwat, satgas ); + + ecl_rft_node_append_cell( rft_node, cell ); + } + + rft ecl_node( rft_node ); + ecl_rft_node_fwrite( ecl_node.get(), fortio.get(), units.getEclType() ); + } + + fortio.close(); +} + +inline std::string uppercase( std::string x ) { + std::transform( x.begin(), x.end(), x.begin(), + []( char c ) { return std::toupper( c ); } ); + + return x; +} + +} + +class EclipseIO::Impl { + public: + Impl( const EclipseState& es, EclipseGrid grid ); + void writeINITFile( const data::Solution& simProps, const NNC& nnc) const; + void writeEGRIDFile( const NNC& nnc ) const; + + EclipseGrid grid; + const EclipseState& es; + std::string outputDir; + std::string baseName; + out::Summary summary; + RFT rft; + bool output_enabled; +}; + +EclipseIO::Impl::Impl( const EclipseState& eclipseState, + EclipseGrid grid_) + : es( eclipseState ) + , grid( std::move( grid_ ) ) + , outputDir( eclipseState.getIOConfig().getOutputDir() ) + , baseName( uppercase( eclipseState.getIOConfig().getBaseName() ) ) + , summary( eclipseState, eclipseState.getSummaryConfig() , grid ) + , rft( outputDir.c_str(), baseName.c_str(), es.getIOConfig().getFMTOUT() ) + , output_enabled( eclipseState.getIOConfig().getOutputEnabled() ) +{} + + + +void EclipseIO::Impl::writeINITFile( const data::Solution& simProps, const NNC& nnc) const { + const auto& units = this->es.getUnits(); + const IOConfig& ioConfig = this->es.cfg().io(); + + std::string initFile( ERT::EclFilename( this->outputDir, + this->baseName, + ECL_INIT_FILE, + ioConfig.getFMTOUT() )); + + ERT::FortIO fortio( initFile, + std::ios_base::out, + ioConfig.getFMTOUT(), + ECL_ENDIAN_FLIP ); + + + // Write INIT header. Observe that the PORV vector is treated + // specially; that is because for this particulat vector we write + // a total of nx*ny*nz values, where the PORV vector has been + // explicitly set to zero for inactive cells. The convention is + // that the active/inactive cell mapping can be inferred by + // reading the PORV vector. + { + + const auto& opm_data = this->es.get3DProperties().getDoubleGridProperty("PORV").getData(); + auto ecl_data = opm_data; + + for (size_t global_index = 0; global_index < opm_data.size(); global_index++) + if (!this->grid.cellActive( global_index )) + ecl_data[global_index] = 0; + + + ecl_init_file_fwrite_header( fortio.get(), + this->grid.c_ptr(), + NULL, + units.getEclType(), + this->es.runspec( ).eclPhaseMask( ), + this->es.getSchedule().posixStartTime( )); + + units.from_si( UnitSystem::measure::volume, ecl_data ); + writeKeyword( fortio, "PORV" , ecl_data ); + } + + // Writing quantities which are calculated by the grid to the INIT file. + ecl_grid_fwrite_depth( this->grid.c_ptr() , fortio.get() , units.getEclType( ) ); + ecl_grid_fwrite_dims( this->grid.c_ptr() , fortio.get() , units.getEclType( ) ); + + // Write properties from the input deck. + { + const auto& properties = this->es.get3DProperties(); + using double_kw = std::pair; + std::vector doubleKeywords = {{"PORO" , UnitSystem::measure::identity }, + {"PERMX" , UnitSystem::measure::permeability }, + {"PERMY" , UnitSystem::measure::permeability }, + {"PERMZ" , UnitSystem::measure::permeability }}; + + for (const auto& kw_pair : doubleKeywords) { + const auto& opm_property = properties.getDoubleGridProperty(kw_pair.first); + auto ecl_data = opm_property.compressedCopy( this->grid ); + + units.from_si( kw_pair.second, ecl_data ); + writeKeyword( fortio, kw_pair.first, ecl_data ); + } + } + + + // Write properties which have been initialized by the simulator. + { + for (const auto& prop : simProps) { + auto ecl_data = this->grid.compressedVector( prop.second.data ); + writeKeyword( fortio, prop.first, ecl_data ); + } + } + + // Write tables + { + Tables tables( this->es.getUnits() ); + tables.addPVTO( this->es.getTableManager().getPvtoTables() ); + tables.addPVTG( this->es.getTableManager().getPvtgTables() ); + tables.fwrite( fortio ); + } + + // Write all integer field properties from the input deck. + { + const auto& properties = this->es.get3DProperties().getIntProperties(); + + // It seems that the INIT file should always contain these + // keywords, we therefor call getKeyword() here to invoke the + // autocreation property, and ensure that the keywords exist + // in the properties container. + properties.getKeyword("PVTNUM"); + properties.getKeyword("SATNUM"); + properties.getKeyword("EQLNUM"); + properties.getKeyword("FIPNUM"); + + for (const auto& property : properties) { + auto ecl_data = property.compressedCopy( this->grid ); + writeKeyword( fortio , property.getKeywordName() , ecl_data ); + } + } + + + // Write NNC transmissibilities + { + std::vector tran; + for( const NNCdata& nd : nnc.nncdata() ) + tran.push_back( nd.trans ); + + units.from_si( UnitSystem::measure::transmissibility , tran ); + writeKeyword( fortio, "TRANNNC" , tran ); + } +} + + +void EclipseIO::Impl::writeEGRIDFile( const NNC& nnc ) const { + const auto& ioConfig = this->es.getIOConfig(); + + std::string egridFile( ERT::EclFilename( this->outputDir, + this->baseName, + ECL_EGRID_FILE, + ioConfig.getFMTOUT() )); + + { + int idx = 0; + auto* ecl_grid = const_cast< ecl_grid_type* >( this->grid.c_ptr() ); + for (const NNCdata& n : nnc.nncdata()) + ecl_grid_add_self_nnc( ecl_grid, n.cell1, n.cell2, idx++); + + ecl_grid_fwrite_EGRID2(ecl_grid, egridFile.c_str(), this->es.getDeckUnitSystem().getEclType() ); + } +} + + +void EclipseIO::writeInitial( data::Solution simProps, const NNC& nnc) { + if( !this->impl->output_enabled ) + return; + + { + const auto& es = this->impl->es; + const IOConfig& ioConfig = es.cfg().io(); + + simProps.convertFromSI( es.getUnits() ); + if( ioConfig.getWriteINITFile() ) + this->impl->writeINITFile( simProps , nnc ); + + if( ioConfig.getWriteEGRIDFile( ) ) + this->impl->writeEGRIDFile( nnc ); + } + + this->impl->summary.set_initial( simProps ); +} + + + +// implementation of the writeTimeStep method +void EclipseIO::writeTimeStep(int report_step, + bool isSubstep, + double secs_elapsed, + data::Solution cells, + data::Wells wells) + { + + if( !this->impl->output_enabled ) + return; + + + const auto& es = this->impl->es; + const auto& grid = this->impl->grid; + const auto& units = es.getUnits(); + const auto& ioConfig = es.getIOConfig(); + const auto& restart = es.cfg().restart(); + + + + const auto& schedule = es.getSchedule(); + + /* + This routine can optionally write RFT and/or restart file; to + be certain that the data correctly converted to output units + the conversion is done once here - and not closer to the actual + use-site. + */ + // Write restart file + if(!isSubstep && restart.getWriteRestartFile(report_step)) + { + std::string filename = ERT::EclFilename( this->impl->outputDir, + this->impl->baseName, + ioConfig.getUNIFOUT() ? ECL_UNIFIED_RESTART_FILE : ECL_RESTART_FILE, + report_step, + ioConfig.getFMTOUT() ); + + RestartIO::save( filename , report_step, secs_elapsed, cells, wells, es , grid); + } + + + const auto unit_type = es.getDeckUnitSystem().getType(); + { + std::vector sched_wells = schedule.getWells( report_step ); + const auto rft_active = [report_step] (const Well* w) { return w->getRFTActive( report_step ) || w->getPLTActive( report_step ); }; + if (std::any_of(sched_wells.begin(), sched_wells.end(), rft_active)) { + this->impl->rft.writeTimeStep( sched_wells, + grid, + report_step, + secs_elapsed + schedule.posixStartTime(), + units.from_si( UnitSystem::measure::time, secs_elapsed ), + units, + cells ); + } + } + + if( isSubstep ) return; + + this->impl->summary.add_timestep( report_step, + secs_elapsed, + this->impl->es, + wells , + cells ); + this->impl->summary.write(); + } + + +std::pair< data::Solution, data::Wells > +EclipseIO::loadRestart(const std::map& keys) const { + const auto& es = this->impl->es; + const auto& grid = this->impl->grid; + const InitConfig& initConfig = es.getInitConfig(); + const auto& ioConfig = es.getIOConfig(); + const int report_step = initConfig.getRestartStep(); + const std::string filename = ioConfig.getRestartFileName( initConfig.getRestartRootName(), + report_step, + false ); + + return RestartIO::load( filename , report_step , keys , es, grid ); +} + +EclipseIO::EclipseIO( const EclipseState& es, EclipseGrid grid) + : impl( new Impl( es, std::move( grid ) ) ) +{ + if( !this->impl->output_enabled ) + return; + { + const auto& outputDir = this->impl->outputDir; + + // make sure that the output directory exists, if not try to create it + if ( !util_entry_exists( outputDir.c_str() ) ) { + util_make_path( outputDir.c_str() ); + } + + if( !util_is_directory( outputDir.c_str() ) ) { + throw std::runtime_error( "The path specified as output directory '" + + outputDir + "' is not a directory"); + } + } +} + + +EclipseIO::~EclipseIO() {} + +} // namespace Opm + + + diff --git a/src/opm/output/eclipse/EclipseReader.cpp b/src/opm/output/eclipse/EclipseReader.cpp index 73e30c828..223a07f1d 100644 --- a/src/opm/output/eclipse/EclipseReader.cpp +++ b/src/opm/output/eclipse/EclipseReader.cpp @@ -53,70 +53,29 @@ namespace { } inline data::Solution restoreSOLUTION( ecl_file_type* file, + const std::map& keys, int numcells, const UnitSystem& units ) { - for( const auto* key : { "PRESSURE", "TEMP", "SWAT", "SGAS" } ) { - if( !ecl_file_has_kw( file, key ) ) + data::Solution sol; + for (const auto& pair : keys) { + const std::string& key = pair.first; + UnitSystem::measure dim = pair.second; + if( !ecl_file_has_kw( file, key.c_str() ) ) throw std::runtime_error("Read of restart file: " "File does not contain " - + std::string( key ) + + key + " data" ); - } - const ecl_kw_type * pres = ecl_file_iget_named_kw( file, "PRESSURE", 0 ); - const ecl_kw_type * temp = ecl_file_iget_named_kw( file, "TEMP", 0 ); - const ecl_kw_type * swat = ecl_file_iget_named_kw( file, "SWAT", 0 ); - const ecl_kw_type * sgas = ecl_file_iget_named_kw( file, "SGAS", 0 ); - - for( const auto& kw : { pres, temp, swat, sgas } ) { - if( ecl_kw_get_size(kw) != numcells) + const ecl_kw_type * ecl_kw = ecl_file_iget_named_kw( file , key.c_str() , 0 ); + if( ecl_kw_get_size(ecl_kw) != numcells) throw std::runtime_error("Restart file: Could not restore " - + std::string( ecl_kw_get_header( kw ) ) + + std::string( ecl_kw_get_header( ecl_kw ) ) + ", mismatched number of cells" ); - } - data::Solution sol; - - { - const auto apply_pressure = [&]( double x ) { - return units.to_si( UnitSystem::measure::pressure, x ); - }; - - const auto apply_temperature = [=]( double x ) { - return units.to_si( UnitSystem::measure::temperature, x ); - }; - - std::vector pressure = double_vector( pres ); - std::vector temperature = double_vector( temp ); - - - std::transform( pressure.begin(), pressure.end(), - pressure.begin(), apply_pressure ); - std::transform( temperature.begin(), temperature.end(), - temperature.begin(), apply_temperature ); - - - - sol.insert( "PRESSURE" , UnitSystem::measure::pressure , pressure , data::TargetType::RESTART_SOLUTION ); - sol.insert( "TEMP" , UnitSystem::measure::temperature , temperature , data::TargetType::RESTART_SOLUTION ); - sol.insert( "SWAT" , UnitSystem::measure::identity , double_vector( swat ) , data::TargetType::RESTART_SOLUTION ); - sol.insert( "SGAS" , UnitSystem::measure::identity , double_vector( sgas ) , data::TargetType::RESTART_SOLUTION ); - - /* optional keywords */ - if( ecl_file_has_kw( file, "RS" ) ) { - const ecl_kw_type * rs_kw = ecl_file_iget_named_kw( file , "RS" , 0); - std::vector rs = double_vector( rs_kw ); - units.to_si( UnitSystem::measure::gas_oil_ratio, rs ); - sol.insert( "RS" , UnitSystem::measure::gas_oil_ratio , rs , data::TargetType::RESTART_SOLUTION ); - } - - if( ecl_file_has_kw( file, "RV" ) ) { - const ecl_kw_type * rv_kw = ecl_file_iget_named_kw( file , "RV" , 0); - std::vector rv = double_vector( rv_kw ); - units.to_si( UnitSystem::measure::oil_gas_ratio, rv ); - sol.insert( "RV" , UnitSystem::measure::oil_gas_ratio , rv , data::TargetType::RESTART_SOLUTION ); - } + std::vector data = double_vector( ecl_kw ); + units.to_si( dim , data ); + sol.insert( key, dim, data , data::TargetType::RESTART_SOLUTION ); } return sol; @@ -129,9 +88,15 @@ data::Wells restore_wells( const double* xwel_data, const int* iwel_data, size_t iwel_data_size, int restart_step, - const std::vector< const Well* > sched_wells, - const std::vector< rt >& phases, - const EclipseGrid& grid ) { + const EclipseState& es ) { + + const auto& sched_wells = es.getSchedule().getWells( restart_step ); + const EclipseGrid& grid = es.getInputGrid( ); + std::vector< rt > phases; + const auto& phase = es.runspec().phases(); + if( phase.active( Phase::WATER ) ) phases.push_back( rt::wat ); + if( phase.active( Phase::OIL ) ) phases.push_back( rt::oil ); + if( phase.active( Phase::GAS ) ) phases.push_back( rt::gas ); const auto well_size = [&]( size_t acc, const Well* w ) { @@ -194,7 +159,7 @@ data::Wells restore_wells( const double* xwel_data, /* should take grid as argument because it may be modified from the simulator */ std::pair< data::Solution, data::Wells > -init_from_restart_file( const EclipseState& es, int numcells ) { +load_from_restart_file( const EclipseState& es, const std::map& keys, int numcells ) { const InitConfig& initConfig = es.getInitConfig(); const auto& ioConfig = es.getIOConfig(); @@ -206,14 +171,6 @@ init_from_restart_file( const EclipseState& es, int numcells ) { restart_step, output); const bool unified = ioConfig.getUNIFIN(); - const auto& sched_wells = es.getSchedule().getWells( restart_step ); - - std::vector< rt > phases; - const auto& phase = es.runspec().phases(); - if( phase.active( Phase::WATER ) ) phases.push_back( rt::wat ); - if( phase.active( Phase::OIL ) ) phases.push_back( rt::oil ); - if( phase.active( Phase::GAS ) ) phases.push_back( rt::gas ); - using ft = ERT::ert_unique_ptr< ecl_file_type, ecl_file_close >; ft file( ecl_file_open( filename.c_str(), 0 ) ); @@ -236,13 +193,11 @@ init_from_restart_file( const EclipseState& es, int numcells ) { const auto iwel_size = ecl_kw_get_size( iwel ); return { - restoreSOLUTION( file.get(), numcells, es.getUnits() ), + restoreSOLUTION( file.get(), keys, numcells, es.getUnits() ), restore_wells( xwel_data, xwel_size, iwel_data, iwel_size, restart_step, - sched_wells, - phases, - es.getInputGrid() ) + es ) }; } diff --git a/src/opm/output/eclipse/EclipseWriter.cpp b/src/opm/output/eclipse/EclipseWriter.cpp deleted file mode 100644 index 7f3066d0b..000000000 --- a/src/opm/output/eclipse/EclipseWriter.cpp +++ /dev/null @@ -1,849 +0,0 @@ -/* - Copyright (c) 2016 Statoil ASA - Copyright (c) 2013-2015 Andreas Lauser - Copyright (c) 2013 SINTEF ICT, Applied Mathematics. - Copyright (c) 2013 Uni Research AS - Copyright (c) 2015 IRIS AS - - 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 - -#include "config.h" - -#include "EclipseWriter.hpp" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include // unique_ptr -#include // move - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#define OPM_XWEL "OPM_XWEL" -#define OPM_IWEL "OPM_IWEL" - -// namespace start here since we don't want the ERT headers in it -namespace Opm { -namespace { - - -inline void convertFromSiTo( std::vector< double >& siValues, - const UnitSystem& units, - UnitSystem::measure m ) { - for (size_t curIdx = 0; curIdx < siValues.size(); ++curIdx) { - siValues[curIdx] = units.from_si( m, siValues[ curIdx ] ); - } -} - -inline int to_ert_welltype( const Well& well, size_t timestep ) { - - if( well.isProducer( timestep ) ) return IWEL_PRODUCER; - - switch( well.getInjectionProperties( timestep ).injectorType ) { - case WellInjector::WATER: return IWEL_WATER_INJECTOR; - case WellInjector::GAS: return IWEL_GAS_INJECTOR; - case WellInjector::OIL: return IWEL_OIL_INJECTOR; - default: return IWEL_UNDOCUMENTED_ZERO; - } -} - - -void writeKeyword( ERT::FortIO& fortio , - const std::string& keywordName, - const std::vector &data ) { - ERT::EclKW< int > kw( keywordName, data ); - kw.fwrite( fortio ); -} - -/* - This overload hardcodes the common assumption that properties which - are stored internally as double values in OPM should be stored as - float values in the ECLIPSE formatted binary files. -*/ - -void writeKeyword( ERT::FortIO& fortio , - const std::string& keywordName, - const std::vector &data) { - - ERT::EclKW< float > kw( keywordName, data ); - kw.fwrite( fortio ); - -} - - -/** - * Pointer to memory that holds the name to an Eclipse output file. - */ -class FileName { -public: - FileName(const std::string& outputDir, - const std::string& baseName, - ecl_file_enum type, - int writeStepIdx, - bool formatted ) : - filename( ecl_util_alloc_filename( - outputDir.c_str(), - baseName.c_str(), - type, - formatted, - writeStepIdx ), - std::free ) - {} - - FileName(const std::string& outputDir, - const std::string& baseName, - ecl_file_enum type, - bool formatted ) : - - filename( ecl_util_alloc_filename( - outputDir.c_str(), - baseName.c_str(), - type, - formatted, - 0 ), - std::free ) - {} - - const char* get() const { return this->filename.get(); } - -private: - using fd = std::unique_ptr< char, decltype( std::free )* >; - fd filename; -}; - -using restart_file = ERT::ert_unique_ptr< ecl_rst_file_type, ecl_rst_file_close >; - -restart_file open_rst( const char* filename, - bool first_restart, - bool unifout, - int report_step ) { - - if( !unifout ) - return restart_file{ ecl_rst_file_open_write( filename ) }; - - if( first_restart ) - return restart_file{ ecl_rst_file_open_write_seek( filename, report_step ) }; - - return restart_file{ ecl_rst_file_open_append( filename ) }; -} - -class Restart { -public: - static const int NIWELZ = 11; //Number of data elements per well in IWEL array in restart file - static const int NZWELZ = 3; //Number of 8-character words per well in ZWEL array restart file - static const int NICONZ = 15; //Number of data elements per completion in ICON array restart file - - /** - * The constants NIWELZ and NZWELZ referes to the number of elements per - * well that we write to the IWEL and ZWEL eclipse restart file data - * arrays. The constant NICONZ refers to the number of elements per - * completion in the eclipse restart file ICON data array.These numbers are - * written to the INTEHEAD header. - - * The elements are added in the method addRestartFileIwelData(...) and and - * addRestartFileIconData(...), respectively. We write as many elements - * that we need to be able to view the restart file in Resinsight. The - * restart file will not be possible to restart from with Eclipse, we write - * to little information to be able to do this. - * - * Observe that all of these values are our "current-best-guess" for how - * many numbers are needed; there might very well be third party - * applications out there which have a hard expectation for these values. - */ - - - Restart(const std::string& outputDir, - const std::string& baseName, - int writeStepIdx, - const IOConfig& ioConfig, - bool first_restart ) : - filename( outputDir, - baseName, - ioConfig.getUNIFOUT() ? ECL_UNIFIED_RESTART_FILE : ECL_RESTART_FILE, - writeStepIdx, - ioConfig.getFMTOUT() ), - rst_file( open_rst( filename.get(), - first_restart, - ioConfig.getUNIFOUT(), - writeStepIdx ) ) - {} - - template< typename T > - void add_kw( ERT::EclKW< T >&& kw ) { - ecl_rst_file_add_kw( this->rst_file.get(), kw.get() ); - } - - void addRestartFileIwelData( std::vector& data, - size_t step, - const Well& well, - size_t offset ) const { - - const auto& completions = well.getCompletions( step ); - - data[ offset + IWEL_HEADI_INDEX ] = well.getHeadI( step ) + 1; - data[ offset + IWEL_HEADJ_INDEX ] = well.getHeadJ( step ) + 1; - data[ offset + IWEL_CONNECTIONS_INDEX ] = completions.size(); - data[ offset + IWEL_GROUP_INDEX ] = 1; - - data[ offset + IWEL_TYPE_INDEX ] = to_ert_welltype( well, step ); - data[ offset + IWEL_STATUS_INDEX ] = - well.getStatus( step ) == WellCommon::OPEN ? 1 : 0; - } - - void addRestartFileIconData( std::vector< int >& data, - const CompletionSet& completions, - size_t wellICONOffset ) const { - - for( size_t i = 0; i < completions.size(); ++i ) { - const auto& completion = completions.get( i ); - size_t offset = wellICONOffset + i * Restart::NICONZ; - data[ offset + ICON_IC_INDEX ] = 1; - - data[ offset + ICON_I_INDEX ] = completion.getI() + 1; - data[ offset + ICON_J_INDEX ] = completion.getJ() + 1; - data[ offset + ICON_K_INDEX ] = completion.getK() + 1; - - const auto open = WellCompletion::StateEnum::OPEN; - data[ offset + ICON_STATUS_INDEX ] = completion.getState() == open - ? 1 - : 0; - - data[ offset + ICON_DIRECTION_INDEX ] = completion.getDirection(); - } - } - - void writeHeader( int stepIdx, ecl_rsthead_type* rsthead_data ) { - ecl_util_set_date_values( rsthead_data->sim_time, - &rsthead_data->day, - &rsthead_data->month, - &rsthead_data->year ); - ecl_rst_file_fwrite_header( this->rst_file.get(), stepIdx, rsthead_data ); - - } - - ecl_rst_file_type* ertHandle() { return this->rst_file.get(); } - const ecl_rst_file_type* ertHandle() const { return this->rst_file.get(); } - -private: - FileName filename; - restart_file rst_file; -}; - -/** - * The Solution class wraps the actions that must be done to the restart file while - * writing solution variables; it is not a handle on its own. - */ -class Solution { -public: - Solution( Restart& res ) : restart( res ) { - ecl_rst_file_start_solution( res.ertHandle() ); - } - - template< typename T > - void add( ERT::EclKW< T >&& kw ) { - ecl_rst_file_add_kw( this->restart.ertHandle(), kw.get() ); - } - - template - void addFromCells(const data::Solution& solution) { - for (const auto& elm: solution) { - if (elm.second.target == data::TargetType::RESTART_SOLUTION) - this->add( ERT::EclKW(elm.first, elm.second.data )); - } - } - - ecl_rst_file_type* ertHandle() { return this->restart.ertHandle(); } - - ~Solution() { ecl_rst_file_end_solution( this->restart.ertHandle() ); } - -private: - Restart& restart; -}; - -/// Convert OPM phase usage to ERT bitmask -inline int ertPhaseMask( const Phases& phase ) { - return ( phase.active( Phase::WATER ) ? ECL_WATER_PHASE : 0 ) - | ( phase.active( Phase::OIL ) ? ECL_OIL_PHASE : 0 ) - | ( phase.active( Phase::GAS ) ? ECL_GAS_PHASE : 0 ); -} - -class RFT { - public: - RFT( const char* output_dir, - const char* basename, - bool format ); - - void writeTimeStep( std::vector< const Well* >, - const EclipseGrid& grid, - int report_step, - time_t current_time, - double days, - ert_ecl_unit_enum, - const data::Solution& cells); - private: - ERT::FortIO fortio; -}; - - -RFT::RFT( const char* output_dir, - const char* basename, - bool format ) : - fortio(FileName( output_dir, basename, ECL_RFT_FILE, format ).get(),std::ios_base::out) -{} - -inline ert_ecl_unit_enum to_ert_unit( UnitSystem::UnitType t ) { - using ut = UnitSystem::UnitType; - switch ( t ) { - case ut::UNIT_TYPE_METRIC: return ECL_METRIC_UNITS; - case ut::UNIT_TYPE_FIELD: return ECL_FIELD_UNITS; - case ut::UNIT_TYPE_LAB: return ECL_LAB_UNITS; - } - - throw std::invalid_argument("unhandled enum value"); -} - -void RFT::writeTimeStep( std::vector< const Well* > wells, - const EclipseGrid& grid, - int report_step, - time_t current_time, - double days, - ert_ecl_unit_enum unitsystem, - const data::Solution& cells) { - - using rft = ERT::ert_unique_ptr< ecl_rft_node_type, ecl_rft_node_free >; - const std::vector& pressure = cells.data("PRESSURE"); - const std::vector& swat = cells.data("SWAT"); - const std::vector& sgas = cells.has("SGAS") ? cells.data("SGAS") : std::vector( pressure.size() , 0 ); - - for( const auto& well : wells ) { - if( !( well->getRFTActive( report_step ) - || well->getPLTActive( report_step ) ) ) - continue; - - auto* rft_node = ecl_rft_node_alloc_new( well->name().c_str(), "RFT", - current_time, days ); - - for( const auto& completion : well->getCompletions( report_step ) ) { - const size_t i = size_t( completion.getI() ); - const size_t j = size_t( completion.getJ() ); - const size_t k = size_t( completion.getK() ); - - if( !grid.cellActive( i, j, k ) ) continue; - - const auto index = grid.activeIndex( i, j, k ); - const double depth = grid.getCellDepth( i, j, k ); - const double press = pressure[ index ]; - const double satwat = swat[ index ]; - const double satgas = sgas[ index ]; - - auto* cell = ecl_rft_cell_alloc_RFT( - i, j, k, depth, press, satwat, satgas ); - - ecl_rft_node_append_cell( rft_node, cell ); - } - - rft ecl_node( rft_node ); - ecl_rft_node_fwrite( ecl_node.get(), this->fortio.get(), unitsystem ); - } -} - -inline std::string uppercase( std::string x ) { - std::transform( x.begin(), x.end(), x.begin(), - []( char c ) { return std::toupper( c ); } ); - - return x; -} - -} - -class EclipseWriter::Impl { - public: - Impl( const EclipseState& es, EclipseGrid grid ); - void writeINITFile( const data::Solution& simProps, const NNC& nnc) const; - void writeEGRIDFile( const NNC& nnc ) const; - - const EclipseState& es; - EclipseGrid grid; - out::RegionCache regionCache; - std::string outputDir; - std::string baseName; - out::Summary summary; - RFT rft; - time_t sim_start_time; - std::array< int, 3 > cartesianSize; - bool output_enabled; - int ert_phase_mask; - bool first_restart = true; -}; - -EclipseWriter::Impl::Impl( const EclipseState& eclipseState, - EclipseGrid grid_) - : es( eclipseState ) - , grid( std::move( grid_ ) ) - , regionCache( es , grid ) - , outputDir( eclipseState.getIOConfig().getOutputDir() ) - , baseName( uppercase( eclipseState.getIOConfig().getBaseName() ) ) - , summary( eclipseState, eclipseState.getSummaryConfig() , grid ) - , rft( outputDir.c_str(), baseName.c_str(), es.getIOConfig().getFMTOUT() ) - , sim_start_time( es.getSchedule().posixStartTime() ) - , output_enabled( eclipseState.getIOConfig().getOutputEnabled() ) - , ert_phase_mask( ertPhaseMask( eclipseState.runspec().phases() ) ) -{} - - - -void EclipseWriter::Impl::writeINITFile( const data::Solution& simProps, const NNC& nnc) const { - const auto& units = this->es.getUnits(); - const IOConfig& ioConfig = this->es.cfg().io(); - - FileName initFile( this->outputDir, - this->baseName, - ECL_INIT_FILE, - ioConfig.getFMTOUT() ); - - ERT::FortIO fortio( initFile.get() , - std::ios_base::out, - ioConfig.getFMTOUT(), - ECL_ENDIAN_FLIP ); - - - // Write INIT header. Observe that the PORV vector is treated - // specially; that is because for this particulat vector we write - // a total of nx*ny*nz values, where the PORV vector has been - // explicitly set to zero for inactive cells. The convention is - // that the active/inactive cell mapping can be inferred by - // reading the PORV vector. - { - - const auto& opm_data = this->es.get3DProperties().getDoubleGridProperty("PORV").getData(); - auto ecl_data = opm_data; - - for (size_t global_index = 0; global_index < opm_data.size(); global_index++) - if (!this->grid.cellActive( global_index )) - ecl_data[global_index] = 0; - - - ecl_init_file_fwrite_header( fortio.get(), - this->grid.c_ptr(), - NULL, - to_ert_unit( units.getType( ) ), - this->ert_phase_mask, - this->sim_start_time); - - convertFromSiTo( ecl_data, - units, - UnitSystem::measure::volume); - - writeKeyword( fortio, "PORV" , ecl_data ); - } - - // Writing quantities which are calculated by the grid to the INIT file. - ecl_grid_fwrite_depth( this->grid.c_ptr() , fortio.get() , to_ert_unit( units.getType( )) ); - ecl_grid_fwrite_dims( this->grid.c_ptr() , fortio.get() , to_ert_unit( units.getType( )) ); - - // Write properties from the input deck. - { - const auto& properties = this->es.get3DProperties(); - using double_kw = std::pair; - std::vector doubleKeywords = {{"PORO" , UnitSystem::measure::identity }, - {"PERMX" , UnitSystem::measure::permeability }, - {"PERMY" , UnitSystem::measure::permeability }, - {"PERMZ" , UnitSystem::measure::permeability }}; - - for (const auto& kw_pair : doubleKeywords) { - const auto& opm_property = properties.getDoubleGridProperty(kw_pair.first); - auto ecl_data = opm_property.compressedCopy( this->grid ); - - convertFromSiTo( ecl_data, - units, - kw_pair.second ); - - writeKeyword( fortio, kw_pair.first, ecl_data ); - } - } - - - // Write properties which have been initialized by the simulator. - { - for (const auto& prop : simProps) { - auto ecl_data = this->grid.compressedVector( prop.second.data ); - writeKeyword( fortio, prop.first, ecl_data ); - } - } - - // Write tables - { - Tables tables( this->es.getUnits() ); - tables.addPVTO( this->es.getTableManager().getPvtoTables() ); - tables.addPVTG( this->es.getTableManager().getPvtgTables() ); - tables.fwrite( fortio ); - } - - // Write all integer field properties from the input deck. - { - const auto& properties = this->es.get3DProperties().getIntProperties(); - - // It seems that the INIT file should always contain these - // keywords, we therefor call getKeyword() here to invoke the - // autocreation property, and ensure that the keywords exist - // in the properties container. - properties.getKeyword("PVTNUM"); - properties.getKeyword("SATNUM"); - properties.getKeyword("EQLNUM"); - properties.getKeyword("FIPNUM"); - - for (const auto& property : properties) { - auto ecl_data = property.compressedCopy( this->grid ); - writeKeyword( fortio , property.getKeywordName() , ecl_data ); - } - } - - - // Write NNC transmissibilities - { - std::vector tran; - for( const NNCdata& nd : nnc.nncdata() ) - tran.push_back( nd.trans ); - - convertFromSiTo( tran, units, UnitSystem::measure::transmissibility ); - writeKeyword( fortio, "TRANNNC" , tran ); - } -} - - -void EclipseWriter::Impl::writeEGRIDFile( const NNC& nnc ) const { - const auto& ioConfig = this->es.getIOConfig(); - - FileName egridFile( this->outputDir, - this->baseName, - ECL_EGRID_FILE, - ioConfig.getFMTOUT() ); - - { - int idx = 0; - auto* ecl_grid = const_cast< ecl_grid_type* >( this->grid.c_ptr() ); - for (const NNCdata& n : nnc.nncdata()) - ecl_grid_add_self_nnc( ecl_grid, n.cell1, n.cell2, idx++); - - ecl_grid_fwrite_EGRID2(ecl_grid, egridFile.get(), to_ert_unit( this->es.getDeckUnitSystem().getType())); - } -} - - -void EclipseWriter::writeInitial( data::Solution simProps, const NNC& nnc) { - if( !this->impl->output_enabled ) - return; - - { - const auto& es = this->impl->es; - const IOConfig& ioConfig = es.cfg().io(); - - simProps.convertFromSI( es.getUnits() ); - if( ioConfig.getWriteINITFile() ) - this->impl->writeINITFile( simProps , nnc ); - - if( ioConfig.getWriteEGRIDFile( ) ) - this->impl->writeEGRIDFile( nnc ); - } - - this->impl->summary.set_initial( simProps ); -} - -std::vector< double > serialize_XWEL( const data::Wells& wells, - int report_step, - const std::vector< const Well* > sched_wells, - const Phases& phase_spec, - const EclipseGrid& grid ) { - - using rt = data::Rates::opt; - - std::vector< rt > phases; - if( phase_spec.active( Phase::WATER ) ) phases.push_back( rt::wat ); - if( phase_spec.active( Phase::OIL ) ) phases.push_back( rt::oil ); - if( phase_spec.active( Phase::GAS ) ) phases.push_back( rt::gas ); - - std::vector< double > xwel; - for( const auto* sched_well : sched_wells ) { - - if( wells.count( sched_well->name() ) == 0 ) { - const auto elems = (sched_well->getCompletions( report_step ).size() - * (phases.size() + data::Completion::restart_size)) - + 2 /* bhp, temperature */ - + phases.size(); - - // write zeros if no well data is provided - xwel.insert( xwel.end(), elems, 0.0 ); - continue; - } - - const auto& well = wells.at( sched_well->name() ); - - xwel.push_back( well.bhp ); - xwel.push_back( well.temperature ); - for( auto phase : phases ) - xwel.push_back( well.rates.get( phase ) ); - - for( const auto& sc : sched_well->getCompletions( report_step ) ) { - const auto i = sc.getI(), j = sc.getJ(), k = sc.getK(); - - const auto rs_size = phases.size() + data::Completion::restart_size; - if( !grid.cellActive( i, j, k ) || sc.getState() == WellCompletion::SHUT ) { - xwel.insert( xwel.end(), rs_size, 0.0 ); - continue; - } - - const auto active_index = grid.activeIndex( i, j, k ); - const auto at_index = [=]( const data::Completion& c ) { - return c.index == active_index; - }; - const auto& completion = std::find_if( well.completions.begin(), - well.completions.end(), - at_index ); - - if( completion == well.completions.end() ) { - xwel.insert( xwel.end(), rs_size, 0.0 ); - continue; - } - - xwel.push_back( completion->pressure ); - xwel.push_back( completion->reservoir_rate ); - for( auto phase : phases ) - xwel.push_back( completion->rates.get( phase ) ); - } - } - - return xwel; -}; - -std::vector< int > serialize_IWEL( const data::Wells& wells, - const std::vector< const Well* > sched_wells ) { - - const auto getctrl = [&]( const Well* w ) { - const auto itr = wells.find( w->name() ); - return itr == wells.end() ? 0 : itr->second.control; - }; - - std::vector< int > iwel( sched_wells.size(), 0.0 ); - std::transform( sched_wells.begin(), sched_wells.end(), iwel.begin(), getctrl ); - - return iwel; -} - -// implementation of the writeTimeStep method -void EclipseWriter::writeTimeStep(int report_step, - bool isSubstep, - double secs_elapsed, - data::Solution cells, - data::Wells wells, - bool write_float) -{ - - if( !this->impl->output_enabled ) - return; - - - time_t current_posix_time = this->impl->sim_start_time + secs_elapsed; - const auto& es = this->impl->es; - const auto& grid = this->impl->grid; - const auto& units = es.getUnits(); - const data::Solution cells_si = cells; - const auto& restart = es.cfg().restart(); - - - const auto days = units.from_si( UnitSystem::measure::time, secs_elapsed ); - const auto& schedule = es.getSchedule(); - - /* - This routine can optionally write RFT and/or restart file; to - be certain that the data correctly converted to output units - the conversion is done once here - and not closer to the actual - use-site. - */ - cells.convertFromSI( units ); - - // Write restart file - if(!isSubstep && restart.getWriteRestartFile(report_step)) - { - const size_t ncwmax = schedule.getMaxNumCompletionsForWells(report_step); - const size_t numWells = schedule.numWells(report_step); - auto wells_ptr = schedule.getWells(report_step); - - std::vector zwell_data( numWells * Restart::NZWELZ , ""); - std::vector iwell_data( numWells * Restart::NIWELZ , 0 ); - std::vector icon_data( numWells * ncwmax * Restart::NICONZ , 0 ); - - Restart restartHandle( this->impl->outputDir, - this->impl->baseName, - report_step, - es.cfg().io(), - this->impl->first_restart ); - - this->impl->first_restart = false; - - for (size_t iwell = 0; iwell < wells_ptr.size(); ++iwell) { - const auto& well = *wells_ptr[iwell]; - { - size_t wellIwelOffset = Restart::NIWELZ * iwell; - restartHandle.addRestartFileIwelData(iwell_data, report_step, well , wellIwelOffset); - } - { - size_t wellIconOffset = ncwmax * Restart::NICONZ * iwell; - restartHandle.addRestartFileIconData(icon_data, well.getCompletions( report_step ), wellIconOffset); - } - zwell_data[ iwell * Restart::NZWELZ ] = well.name().c_str(); - } - - - { - ecl_rsthead_type rsthead_data = {}; - rsthead_data.sim_time = current_posix_time; - rsthead_data.nactive = grid.getNumActive(), - rsthead_data.nx = grid.getNX(); - rsthead_data.ny = grid.getNY(); - rsthead_data.nz = grid.getNZ(); - rsthead_data.nwells = numWells; - rsthead_data.niwelz = Restart::NIWELZ; - rsthead_data.nzwelz = Restart::NZWELZ; - rsthead_data.niconz = Restart::NICONZ; - rsthead_data.ncwmax = ncwmax; - rsthead_data.phase_sum = this->impl->ert_phase_mask; - rsthead_data.sim_days = days; - rsthead_data.unit_system= to_ert_unit( units.getType() ); - - restartHandle.writeHeader( report_step, &rsthead_data); - } - - const auto& phases = es.runspec().phases(); - const auto& sched_wells = schedule.getWells( report_step ); - const auto xwel = serialize_XWEL( wells, report_step, sched_wells, phases, grid ); - const auto iwel = serialize_IWEL( wells, sched_wells ); - - restartHandle.add_kw( ERT::EclKW< int >(IWEL_KW, iwell_data) ); - restartHandle.add_kw( ERT::EclKW< const char* >(ZWEL_KW, zwell_data ) ); - restartHandle.add_kw( ERT::EclKW< double >( OPM_XWEL, xwel ) ); - restartHandle.add_kw( ERT::EclKW< int >( OPM_IWEL, iwel ) ); - restartHandle.add_kw( ERT::EclKW< int >( ICON_KW, icon_data ) ); - - - /* - The code in the block below is in a state of flux, and not - very well tested. In particular there have been issues with - the mapping between active an inactive cells in the past - - there might be more those. - - Currently the code hard-codes the following assumptions: - - 1. All the cells[ ] vectors have size equal to the number - of active cells, and they are already in eclipse - ordering. - - 2. No unit transformatio applied to the saturation and - RS/RV vectors. - - */ - { - Solution sol(restartHandle); //Type solution: confusing - - if (write_float) - sol.addFromCells( cells ); - else - sol.addFromCells( cells ); - // MIssing a ENDSOL output. - { - for (const auto& prop : cells) { - if (prop.second.target == data::TargetType::RESTART_AUXILLARY) { - auto ecl_data = grid.compressedVector( prop.second.data ); - - if (write_float) - sol.add( ERT::EclKW(prop.first, ecl_data)); - else - sol.add( ERT::EclKW(prop.first, ecl_data)); - } - } - } - } - } - - const auto unit_type = es.getDeckUnitSystem().getType(); - this->impl->rft.writeTimeStep( schedule.getWells( report_step ), - grid, - report_step, - current_posix_time, - days, - to_ert_unit( unit_type ), - cells ); - - if( isSubstep ) return; - - this->impl->summary.add_timestep( report_step, - secs_elapsed, - this->impl->es, - this->impl->regionCache, - wells , - cells_si ); - this->impl->summary.write(); - } - - - -EclipseWriter::EclipseWriter( const EclipseState& es, EclipseGrid grid) - : impl( new Impl( es, std::move( grid ) ) ) -{ - if( !this->impl->output_enabled ) - return; - { - const auto& outputDir = this->impl->outputDir; - - // make sure that the output directory exists, if not try to create it - if ( !util_entry_exists( outputDir.c_str() ) ) { - util_make_path( outputDir.c_str() ); - } - - if( !util_is_directory( outputDir.c_str() ) ) { - throw std::runtime_error( "The path specified as output directory '" - + outputDir + "' is not a directory"); - } - } -} - -EclipseWriter::~EclipseWriter() {} - -} // namespace Opm diff --git a/src/opm/output/eclipse/RestartIO.cpp b/src/opm/output/eclipse/RestartIO.cpp new file mode 100644 index 000000000..2e38d106e --- /dev/null +++ b/src/opm/output/eclipse/RestartIO.cpp @@ -0,0 +1,514 @@ +/* + Copyright (c) 2016 Statoil ASA + Copyright (c) 2013-2015 Andreas Lauser + Copyright (c) 2013 SINTEF ICT, Applied Mathematics. + Copyright (c) 2013 Uni Research AS + Copyright (c) 2015 IRIS AS + + 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 +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define OPM_XWEL "OPM_XWEL" +#define OPM_IWEL "OPM_IWEL" + +namespace Opm { +namespace RestartIO { + + + +namespace { + + static const int NIWELZ = 11; //Number of data elements per well in IWEL array in restart file + static const int NZWELZ = 3; //Number of 8-character words per well in ZWEL array restart file + static const int NICONZ = 15; //Number of data elements per completion in ICON array restart file + + /** + * The constants NIWELZ and NZWELZ referes to the number of + * elements per well that we write to the IWEL and ZWEL eclipse + * restart file data arrays. The constant NICONZ refers to the + * number of elements per completion in the eclipse restart file + * ICON data array.These numbers are written to the INTEHEAD + * header. + * + * Observe that all of these values are our "current-best-guess" + * for how many numbers are needed; there might very well be third + * party applications out there which have a hard expectation for + * these values. + */ + + + inline int to_ert_welltype( const Well& well, size_t timestep ) { + if( well.isProducer( timestep ) ) return IWEL_PRODUCER; + + switch( well.getInjectionProperties( timestep ).injectorType ) { + case WellInjector::WATER: + return IWEL_WATER_INJECTOR; + case WellInjector::GAS: + return IWEL_GAS_INJECTOR; + case WellInjector::OIL: + return IWEL_OIL_INJECTOR; + default: + return IWEL_UNDOCUMENTED_ZERO; + } + } + + std::vector double_vector( const ecl_kw_type * ecl_kw ) { + size_t size = ecl_kw_get_size( ecl_kw ); + + if (ecl_kw_get_type( ecl_kw ) == ECL_DOUBLE_TYPE ) { + const double * ecl_data = ecl_kw_get_double_ptr( ecl_kw ); + return { ecl_data , ecl_data + size }; + } else { + const float * ecl_data = ecl_kw_get_float_ptr( ecl_kw ); + return { ecl_data , ecl_data + size }; + } + + } + + + inline data::Solution restoreSOLUTION( ecl_file_view_type* file_view, + const std::map& keys, + const UnitSystem& units, + int numcells) { + + data::Solution sol; + for (const auto& pair : keys) { + const std::string& key = pair.first; + UnitSystem::measure dim = pair.second; + if( !ecl_file_view_has_kw( file_view, key.c_str() ) ) + throw std::runtime_error("Read of restart file: " + "File does not contain " + + key + + " data" ); + + const ecl_kw_type * ecl_kw = ecl_file_view_iget_named_kw( file_view , key.c_str() , 0 ); + if( ecl_kw_get_size(ecl_kw) != numcells) + throw std::runtime_error("Restart file: Could not restore " + + std::string( ecl_kw_get_header( ecl_kw ) ) + + ", mismatched number of cells" ); + + std::vector data = double_vector( ecl_kw ); + units.to_si( dim , data ); + + sol.insert( key, dim, data , data::TargetType::RESTART_SOLUTION ); + } + + return sol; + } + + +using rt = data::Rates::opt; +data::Wells restore_wells( const ecl_kw_type * opm_xwel, + const ecl_kw_type * opm_iwel, + const ecl_kw_type * zwel, + int restart_step, + const EclipseState& es ) { + + const auto& sched_wells = es.getSchedule().getWells( restart_step ); + const EclipseGrid& grid = es.getInputGrid( ); + std::vector< rt > phases; + const auto& phase = es.runspec().phases(); + if( phase.active( Phase::WATER ) ) phases.push_back( rt::wat ); + if( phase.active( Phase::OIL ) ) phases.push_back( rt::oil ); + if( phase.active( Phase::GAS ) ) phases.push_back( rt::gas ); + + + const auto well_size = [&]( size_t acc, const Well* w ) { + return acc + + 2 + phases.size() + + ( w->getCompletions( restart_step ).size() + * (phases.size() + data::Completion::restart_size) ); + }; + + const auto expected_xwel_size = std::accumulate( sched_wells.begin(), + sched_wells.end(), + size_t( 0 ), + well_size ); + + if( ecl_kw_get_size( opm_xwel ) != expected_xwel_size ) { + throw std::runtime_error( + "Mismatch between OPM_XWEL and deck; " + "OPM_XWEL size was " + std::to_string( ecl_kw_get_size( opm_xwel ) ) + + ", expected " + std::to_string( expected_xwel_size ) ); + } + + if( ecl_kw_get_size( opm_iwel ) != sched_wells.size() ) + throw std::runtime_error( + "Mismatch between OPM_IWEL and deck; " + "OPM_IWEL size was " + std::to_string( ecl_kw_get_size( opm_iwel ) ) + + ", expected " + std::to_string( sched_wells.size() ) ); + + data::Wells wells; + const double * opm_xwel_data = ecl_kw_get_double_ptr( opm_xwel ); + const int * opm_iwel_data = ecl_kw_get_int_ptr( opm_iwel ); + for( const auto* sched_well : sched_wells ) { + data::Well& well = wells[ sched_well->name() ]; + + well.bhp = *opm_xwel_data++; + well.temperature = *opm_xwel_data++; + well.control = *opm_iwel_data++; + + for( auto phase : phases ) + well.rates.set( phase, *opm_xwel_data++ ); + + for( const auto& sc : sched_well->getCompletions( restart_step ) ) { + const auto i = sc.getI(), j = sc.getJ(), k = sc.getK(); + if( !grid.cellActive( i, j, k ) || sc.getState() == WellCompletion::SHUT ) { + opm_xwel_data += data::Completion::restart_size + phases.size(); + continue; + } + + const auto active_index = grid.activeIndex( i, j, k ); + + well.completions.emplace_back(); + auto& completion = well.completions.back(); + completion.index = active_index; + completion.pressure = *opm_xwel_data++; + completion.reservoir_rate = *opm_xwel_data++; + for( auto phase : phases ) + completion.rates.set( phase, *opm_xwel_data++ ); + } + } + + return wells; +} +} + +/* should take grid as argument because it may be modified from the simulator */ +std::pair< data::Solution, data::Wells > load( const std::string& filename, + int report_step, + const std::map& keys, + const EclipseState& es, + const EclipseGrid& grid) { + + const bool unified = ( ERT::EclFiletype( filename ) == ECL_UNIFIED_RESTART_FILE ); + ERT::ert_unique_ptr< ecl_file_type, ecl_file_close > file(ecl_file_open( filename.c_str(), 0 )); + ecl_file_view_type * file_view; + + if( !file ) + throw std::runtime_error( "Restart file " + filename + " not found!" ); + + if( unified ) { + file_view = ecl_file_get_restart_view( file.get() , -1 , report_step , -1 , -1 ); + if (!file_view) + throw std::runtime_error( "Restart file " + filename + + " does not contain data for report step " + + std::to_string( report_step ) + "!" ); + } else + file_view = ecl_file_get_global_view( file.get() ); + + const ecl_kw_type * intehead = ecl_file_view_iget_named_kw( file_view , "INTEHEAD", 0 ); + const ecl_kw_type * opm_xwel = ecl_file_view_iget_named_kw( file_view , "OPM_XWEL", 0 ); + const ecl_kw_type * opm_iwel = ecl_file_view_iget_named_kw( file_view, "OPM_IWEL", 0 ); + const ecl_kw_type * zwel = ecl_file_view_iget_named_kw( file_view , "ZWEL" , 0 ); + + UnitSystem units( static_cast(ecl_kw_iget_int( intehead , INTEHEAD_UNIT_INDEX ))); + return { + restoreSOLUTION( file_view, keys, units , grid.getNumActive( )), + restore_wells( opm_xwel, opm_iwel,zwel, + report_step, + es ) + }; +} + + + + + +namespace { + +std::vector serialize_ICON( int report_step, + int ncwmax, + const std::vector& sched_wells) { + + size_t well_offset = 0; + std::vector data( sched_wells.size() * ncwmax * NICONZ , 0 ); + for (const Well* well : sched_wells) { + const auto& completions = well->getCompletions( report_step ); + size_t completion_offset = 0; + for( const auto& completion : completions) { + size_t offset = well_offset + completion_offset; + data[ offset + ICON_IC_INDEX ] = 1; + + data[ offset + ICON_I_INDEX ] = completion.getI() + 1; + data[ offset + ICON_J_INDEX ] = completion.getJ() + 1; + data[ offset + ICON_K_INDEX ] = completion.getK() + 1; + data[ offset + ICON_DIRECTION_INDEX ] = completion.getDirection(); + { + const auto open = WellCompletion::StateEnum::OPEN; + data[ offset + ICON_STATUS_INDEX ] = completion.getState() == open + ? 1 + : 0; + } + + completion_offset += NICONZ; + } + well_offset += ncwmax * NICONZ; + } + return data; +} + +std::vector serialize_IWEL( size_t step, + const std::vector& wells) { + + std::vector data( wells.size() * NIWELZ , 0 ); + size_t offset = 0; + for (const auto well : wells) { + const auto& completions = well->getCompletions( step ); + + data[ offset + IWEL_HEADI_INDEX ] = well->getHeadI( step ) + 1; + data[ offset + IWEL_HEADJ_INDEX ] = well->getHeadJ( step ) + 1; + data[ offset + IWEL_CONNECTIONS_INDEX ] = completions.size(); + data[ offset + IWEL_GROUP_INDEX ] = 1; + + data[ offset + IWEL_TYPE_INDEX ] = to_ert_welltype( *well, step ); + data[ offset + IWEL_STATUS_INDEX ] = + well->getStatus( step ) == WellCommon::OPEN ? 1 : 0; + + offset += NIWELZ; + } + return data; +} + + + + + + + + + + +std::vector< int > serialize_OPM_IWEL( const data::Wells& wells, + const std::vector< const Well* > sched_wells ) { + + const auto getctrl = [&]( const Well* w ) { + const auto itr = wells.find( w->name() ); + return itr == wells.end() ? 0 : itr->second.control; + }; + + std::vector< int > iwel( sched_wells.size(), 0.0 ); + std::transform( sched_wells.begin(), sched_wells.end(), iwel.begin(), getctrl ); + return iwel; +} + +std::vector< double > serialize_OPM_XWEL( const data::Wells& wells, + int report_step, + const std::vector< const Well* > sched_wells, + const Phases& phase_spec, + const EclipseGrid& grid ) { + + using rt = data::Rates::opt; + + std::vector< rt > phases; + if( phase_spec.active( Phase::WATER ) ) phases.push_back( rt::wat ); + if( phase_spec.active( Phase::OIL ) ) phases.push_back( rt::oil ); + if( phase_spec.active( Phase::GAS ) ) phases.push_back( rt::gas ); + + std::vector< double > xwel; + for( const auto* sched_well : sched_wells ) { + + if( wells.count( sched_well->name() ) == 0 ) { + const auto elems = (sched_well->getCompletions( report_step ).size() + * (phases.size() + data::Completion::restart_size)) + + 2 /* bhp, temperature */ + + phases.size(); + + // write zeros if no well data is provided + xwel.insert( xwel.end(), elems, 0.0 ); + continue; + } + + const auto& well = wells.at( sched_well->name() ); + + xwel.push_back( well.bhp ); + xwel.push_back( well.temperature ); + for( auto phase : phases ) + xwel.push_back( well.rates.get( phase ) ); + + for( const auto& sc : sched_well->getCompletions( report_step ) ) { + const auto i = sc.getI(), j = sc.getJ(), k = sc.getK(); + + const auto rs_size = phases.size() + data::Completion::restart_size; + if( !grid.cellActive( i, j, k ) || sc.getState() == WellCompletion::SHUT ) { + xwel.insert( xwel.end(), rs_size, 0.0 ); + continue; + } + + const auto active_index = grid.activeIndex( i, j, k ); + const auto at_index = [=]( const data::Completion& c ) { + return c.index == active_index; + }; + const auto& completion = std::find_if( well.completions.begin(), + well.completions.end(), + at_index ); + + if( completion == well.completions.end() ) { + xwel.insert( xwel.end(), rs_size, 0.0 ); + continue; + } + + xwel.push_back( completion->pressure ); + xwel.push_back( completion->reservoir_rate ); + for( auto phase : phases ) + xwel.push_back( completion->rates.get( phase ) ); + } + } + + return xwel; +}; + + +std::vector serialize_ZWEL( const std::vector& wells) { + std::vector data( wells.size( ) * NZWELZ , ""); + size_t offset = 0; + + for (const auto& well : wells) { + data[ offset ] = well->name().c_str(); + offset += NZWELZ; + } + return data; +} + + + + + +template< typename T > +void write_kw(ecl_rst_file_type * rst_file , ERT::EclKW< T >&& kw) { + ecl_rst_file_add_kw( rst_file, kw.get() ); +} + +void writeHeader(ecl_rst_file_type * rst_file, + int report_step, + time_t posix_time, + double sim_days, + int ert_phase_mask, + const UnitSystem& units, + const Schedule& schedule, + const EclipseGrid& grid) { + + ecl_rsthead_type rsthead_data = {}; + + rsthead_data.sim_time = posix_time; + rsthead_data.nactive = grid.getNumActive(); + rsthead_data.nx = grid.getNX(); + rsthead_data.ny = grid.getNY(); + rsthead_data.nz = grid.getNZ(); + rsthead_data.nwells = schedule.numWells(report_step); + rsthead_data.niwelz = NIWELZ; + rsthead_data.nzwelz = NZWELZ; + rsthead_data.niconz = NICONZ; + rsthead_data.ncwmax = schedule.getMaxNumCompletionsForWells(report_step); + rsthead_data.phase_sum = ert_phase_mask; + rsthead_data.sim_days = sim_days; + rsthead_data.unit_system = units.getEclType( ); + + ecl_util_set_date_values( rsthead_data.sim_time, + &rsthead_data.day, + &rsthead_data.month, + &rsthead_data.year ); + + ecl_rst_file_fwrite_header( rst_file, report_step , &rsthead_data ); +} + +void writeSolution(ecl_rst_file_type* rst_file, const data::Solution& solution) { + ecl_rst_file_start_solution( rst_file ); + for (const auto& elm: solution) { + if (elm.second.target == data::TargetType::RESTART_SOLUTION) + ecl_rst_file_add_kw( rst_file , ERT::EclKW(elm.first, elm.second.data).get()); + } + ecl_rst_file_end_solution( rst_file ); + + for (const auto& elm: solution) { + if (elm.second.target == data::TargetType::RESTART_AUXILLARY) + ecl_rst_file_add_kw( rst_file , ERT::EclKW(elm.first, elm.second.data).get()); + } +} + + +void writeWell(ecl_rst_file_type* rst_file, int report_step, const EclipseState& es , const EclipseGrid& grid, const data::Wells& wells) { + + const auto& schedule = es.getSchedule(); + const auto sched_wells = schedule.getWells(report_step); + const auto& phases = es.runspec().phases(); + const size_t ncwmax = schedule.getMaxNumCompletionsForWells(report_step); + + const auto opm_xwel = serialize_OPM_XWEL( wells, report_step, sched_wells, phases, grid ); + const auto opm_iwel = serialize_OPM_IWEL( wells, sched_wells ); + const auto iwel_data = serialize_IWEL(report_step, sched_wells); + const auto icon_data = serialize_ICON(report_step , ncwmax, sched_wells); + const auto zwel_data = serialize_ZWEL( sched_wells ); + + write_kw( rst_file, ERT::EclKW< int >( IWEL_KW, iwel_data) ); + write_kw( rst_file, ERT::EclKW< const char* >(ZWEL_KW, zwel_data ) ); + write_kw( rst_file, ERT::EclKW< double >( OPM_XWEL, opm_xwel ) ); + write_kw( rst_file, ERT::EclKW< int >( OPM_IWEL, opm_iwel ) ); + write_kw( rst_file, ERT::EclKW< int >( ICON_KW, icon_data ) ); +} +} + + +void save(const std::string& filename, + int report_step, + double seconds_elapsed, + data::Solution cells, + data::Wells wells, + const EclipseState& es, + const EclipseGrid& grid) +{ + const auto& ioConfig = es.getIOConfig(); + int ert_phase_mask = es.runspec().eclPhaseMask( ); + const Schedule& schedule = es.getSchedule(); + const auto& units = es.getUnits(); + time_t posix_time = schedule.posixStartTime() + seconds_elapsed; + const auto sim_time = units.from_si( UnitSystem::measure::time, seconds_elapsed ); + ERT::ert_unique_ptr< ecl_rst_file_type, ecl_rst_file_close > rst_file; + + if (ERT::EclFiletype( filename ) == ECL_UNIFIED_RESTART_FILE) + rst_file.reset( ecl_rst_file_open_write_seek( filename.c_str(), report_step ) ); + else + rst_file.reset( ecl_rst_file_open_write( filename.c_str() ) ); + + cells.convertFromSI( units ); + writeHeader( rst_file.get() , report_step, posix_time , sim_time, ert_phase_mask, units, schedule , grid ); + writeWell( rst_file.get() , report_step, es , grid, wells); + writeSolution( rst_file.get() , cells ); + +} +} +} + diff --git a/src/opm/output/eclipse/Summary.cpp b/src/opm/output/eclipse/Summary.cpp index c4b09cd61..8978652c5 100644 --- a/src/opm/output/eclipse/Summary.cpp +++ b/src/opm/output/eclipse/Summary.cpp @@ -783,6 +783,7 @@ Summary::Summary( const EclipseState& st, const EclipseGrid& grid_arg, const char* basename ) : grid( grid_arg ), + regionCache( st , grid_arg ), ecl_sum( ecl_sum_alloc_writer( basename, @@ -837,7 +838,6 @@ Summary::Summary( const EclipseState& st, void Summary::add_timestep( int report_step, double secs_elapsed, const EclipseState& es, - const RegionCache& regionCache, const data::Wells& wells , const data::Solution& state) { @@ -858,7 +858,7 @@ void Summary::add_timestep( int report_step, num, wells, state, - regionCache, + this->regionCache, this->grid, this->initial_oip } ); diff --git a/tests/test_EclipseWriter.cpp b/tests/test_EclipseIO.cpp similarity index 98% rename from tests/test_EclipseWriter.cpp rename to tests/test_EclipseIO.cpp index f22579c3b..9f6a00969 100644 --- a/tests/test_EclipseWriter.cpp +++ b/tests/test_EclipseIO.cpp @@ -22,10 +22,10 @@ #define BOOST_TEST_DYN_LINK #endif -#define BOOST_TEST_MODULE EclipseWriter +#define BOOST_TEST_MODULE EclipseIO #include -#include +#include #include #include @@ -109,7 +109,7 @@ void compareErtData(const std::vector &src, const std::vector &dst) } void checkEgridFile( const EclipseGrid& eclGrid ) { - // use ERT directly to inspect the EGRID file produced by EclipseWriter + // use ERT directly to inspect the EGRID file produced by EclipseIO auto egridFile = fortio_open_reader("FOO.EGRID", /*isFormated=*/0, ECL_ENDIAN_FLIP); const auto numCells = eclGrid.getNX() * eclGrid.getNY() * eclGrid.getNZ(); @@ -146,7 +146,7 @@ void checkEgridFile( const EclipseGrid& eclGrid ) { } void checkInitFile( const Deck& deck, const data::Solution& simProps) { - // use ERT directly to inspect the INIT file produced by EclipseWriter + // use ERT directly to inspect the INIT file produced by EclipseIO ERT::ert_unique_ptr initFile(ecl_file_open( "FOO.INIT" , 0 )); for (int k=0; k < ecl_file_get_size( initFile.get() ); k++) { @@ -184,7 +184,7 @@ void checkRestartFile( int timeStepIdx ) { for (int i = 1; i <= timeStepIdx; ++i) { auto sol = createBlackoilState( i, 3 * 3 * 3 ); - // use ERT directly to inspect the restart file produced by EclipseWriter + // use ERT directly to inspect the restart file produced by EclipseIO auto rstFile = fortio_open_reader("FOO.UNRST", /*isFormated=*/0, ECL_ENDIAN_FLIP); int curSeqnum = -1; @@ -226,7 +226,7 @@ void checkRestartFile( int timeStepIdx ) { } } -BOOST_AUTO_TEST_CASE(EclipseWriterIntegration) { +BOOST_AUTO_TEST_CASE(EclipseIOIntegration) { const char *deckString = "RUNSPEC\n" "UNIFOUT\n" @@ -276,7 +276,7 @@ BOOST_AUTO_TEST_CASE(EclipseWriterIntegration) { auto& eclGrid = es.getInputGrid(); es.getIOConfig().setBaseName( "FOO" ); - EclipseWriter eclWriter( es, eclGrid ); + EclipseIO eclWriter( es, eclGrid ); using measure = UnitSystem::measure; using TargetType = data::TargetType; diff --git a/tests/test_RFT.cpp b/tests/test_RFT.cpp index c26b0b48d..3d1e0df18 100755 --- a/tests/test_RFT.cpp +++ b/tests/test_RFT.cpp @@ -26,7 +26,7 @@ #include #include -#include +#include #include #include @@ -120,7 +120,7 @@ BOOST_AUTO_TEST_CASE(test_RFT) { const auto& grid = eclipseState.getInputGrid(); const auto numCells = grid.getCartesianSize( ); - EclipseWriter eclipseWriter( eclipseState, grid); + EclipseIO eclipseWriter( eclipseState, grid); time_t start_time = eclipseState.getSchedule().posixStartTime(); /* step time read from deck and hard-coded here */ time_t step_time = ecl_util_make_date(10, 10, 2008 ); diff --git a/tests/test_Restart.cpp b/tests/test_Restart.cpp index 6d357fe79..2b4142cb5 100644 --- a/tests/test_Restart.cpp +++ b/tests/test_Restart.cpp @@ -1,4 +1,3 @@ - /* Copyright 2014 Statoil IT This file is part of the Open Porous Media project (OPM). @@ -22,11 +21,11 @@ #define BOOST_TEST_DYN_LINK #endif -#define BOOST_TEST_MODULE EclipseWriter +#define BOOST_TEST_MODULE EclipseIO #include -#include -#include +#include +#include #include #include #include @@ -45,7 +44,7 @@ #include #include #include -#include +#include using namespace Opm; @@ -266,28 +265,6 @@ bool operator==( const Well& lhs, const Well& rhs ) { } -/* - * forward declarations of internal functions that we want to expose to tests - * but not to users. - */ -std::vector< double > serialize_XWEL( const data::Wells& wells, - int report_step, - const std::vector< const Well* > sched_wells, - const Phases&, - const EclipseGrid& ); - -std::vector< int > serialize_IWEL( const data::Wells& wells, - const std::vector< const Well* > sched_wells ); - -data::Wells restore_wells( const double* xwel_data, - size_t xwel_data_size, - const int* iwel_data, - size_t iwel_data_size, - int restart_step, - const std::vector< const Well* > sched_wells, - const std::vector< data::Rates::opt >& phases, - const EclipseGrid& grid ); -} data::Wells mkWells() { data::Rates r1, r2, rc1, rc2, rc3; @@ -368,17 +345,10 @@ data::Solution mkSolution( int numCells ) { } std::pair< data::Solution, data::Wells > -first_sim(test_work_area_type * test_area) { - - std::string eclipse_data_filename = "FIRST_SIM.DATA"; - test_work_area_copy_file(test_area, eclipse_data_filename.c_str()); - - auto eclipseState = Parser::parse( eclipse_data_filename ); - - const auto& grid = eclipseState.getInputGrid(); +first_sim(const EclipseState& es, EclipseIO& eclWriter) { + const auto& grid = es.getInputGrid(); auto num_cells = grid.getNX() * grid.getNY() * grid.getNZ(); - EclipseWriter eclWriter( eclipseState, grid); auto start_time = ecl_util_make_date( 1, 11, 1979 ); auto first_step = ecl_util_make_date( 10, 10, 2008 ); @@ -393,25 +363,21 @@ first_sim(test_work_area_type * test_area) { return { sol, wells }; } -std::pair< data::Solution, data::Wells > second_sim() { - auto eclipseState = Parser::parseData( input() ); - - const auto& grid = eclipseState.getInputGrid(); - auto num_cells = grid.getNX() * grid.getNY() * grid.getNZ(); - - return init_from_restart_file( eclipseState, num_cells ); +std::pair< data::Solution, data::Wells > second_sim(const EclipseIO& writer, const std::map& keys) { + return writer.loadRestart( keys ); } + void compare( std::pair< data::Solution, data::Wells > fst, - std::pair< data::Solution, data::Wells > snd ) { + std::pair< data::Solution, data::Wells > snd , + const std::map& keys) { - for( auto key : { "PRESSURE", "TEMP", "SWAT", "SGAS", - "RS", "RV" } ) { + for (const auto& pair : keys) { - auto first = fst.first.data( key ).begin(); - auto second = snd.first.data( key ).begin(); + auto first = fst.first.data( pair.first ).begin(); + auto second = snd.first.data( pair.first ).begin(); - for( ; first != fst.first.data( key ).end(); ++first, ++second ) + for( ; first != fst.first.data( pair.first ).end(); ++first, ++second ) BOOST_CHECK_CLOSE( *first, *second, 0.00001 ); } @@ -419,38 +385,23 @@ void compare( std::pair< data::Solution, data::Wells > fst, } BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData) { - test_work_area_type * test_area = test_work_area_alloc("test_Restart"); + std::map keys {{"PRESSURE" , UnitSystem::measure::pressure}, + {"SWAT" , UnitSystem::measure::identity}, + {"SGAS" , UnitSystem::measure::identity}, + {"TEMP" , UnitSystem::measure::temperature}}; + ERT::TestArea testArea("test_Restart"); + testArea.copyFile( "FIRST_SIM.DATA" ); - auto state1 = first_sim(test_area); - auto state2 = second_sim(); - compare(state1, state2); + auto eclipseState = Parser::parse( "FIRST_SIM.DATA" ); + const auto& grid = eclipseState.getInputGrid(); + EclipseIO eclWriter( eclipseState, grid); + auto state1 = first_sim( eclipseState , eclWriter ); + auto state2 = second_sim( eclWriter , keys ); + compare(state1, state2 , keys); - test_work_area_free(test_area); + { + std::map extra_keys {{"SOIL" , UnitSystem::measure::pressure}}; + BOOST_CHECK_THROW( second_sim( eclWriter, extra_keys ) , std::runtime_error ); + } } - -BOOST_AUTO_TEST_CASE(OPM_XWEL) { - auto es = Parser::parseData( input( "XWEL" ) ); - const auto& sched = es.getSchedule(); - const auto& grid = es.getInputGrid(); - const auto& ph = es.runspec().phases(); - - std::vector< data::Rates::opt > phases { - data::Rates::opt::wat, - data::Rates::opt::oil, - data::Rates::opt::gas, - }; - - const auto wells = mkWells(); - const auto& sched_wells = sched.getWells( 1 ); - const auto xwel = serialize_XWEL( wells, 1, sched_wells, ph, grid ); - const auto iwel = serialize_IWEL( wells, sched_wells ); - - const auto restored_wells = restore_wells( xwel.data(), xwel.size(), - iwel.data(), iwel.size(), - 1, - sched.getWells( 1 ), - phases, - grid ); - - BOOST_CHECK_EQUAL( wells, restored_wells ); } diff --git a/tests/test_Summary.cpp b/tests/test_Summary.cpp index a930bceb6..e189bdc15 100644 --- a/tests/test_Summary.cpp +++ b/tests/test_Summary.cpp @@ -199,7 +199,6 @@ struct setup { EclipseState es; SummaryConfig config; const EclipseGrid& grid; - const out::RegionCache regionCache; data::Wells wells; std::string name; ERT::TestArea ta; @@ -214,7 +213,6 @@ struct setup { es( deck, ParseContext() ), config( deck, es, parseContext ), grid( es.getInputGrid() ), - regionCache( es , grid ), wells( result_wells() ), name( fname ), ta( ERT::TestArea("test_summary") ), @@ -239,9 +237,9 @@ BOOST_AUTO_TEST_CASE(well_keywords) { cfg.name = "PATH/CASE"; out::Summary writer( cfg.es, cfg.config, cfg.grid , cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -375,9 +373,9 @@ BOOST_AUTO_TEST_CASE(group_keywords) { setup cfg( "test_Summary_group" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -473,9 +471,9 @@ BOOST_AUTO_TEST_CASE(completion_kewords) { setup cfg( "test_Summary_completion" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -526,9 +524,9 @@ BOOST_AUTO_TEST_CASE(field_keywords) { setup cfg( "test_Summary_field" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -623,9 +621,9 @@ BOOST_AUTO_TEST_CASE(foe_test) { out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); writer.set_initial( sol ); - writer.add_timestep( 1, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 5 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 10 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 1, 2 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 5 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 10 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -640,9 +638,9 @@ BOOST_AUTO_TEST_CASE(report_steps_time) { setup cfg( "test_Summary_report_steps_time" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 1, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 5 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 10 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 1, 2 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 5 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 10 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -662,9 +660,9 @@ BOOST_AUTO_TEST_CASE(skip_unknown_var) { setup cfg( "test_Summary_skip_unknown_var" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 1, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 5 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 10 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 1, 2 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 5 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 10 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -682,9 +680,9 @@ BOOST_AUTO_TEST_CASE(region_vars) { { out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 1, 2 * day, cfg.es, cfg.regionCache, cfg.wells, cfg.solution); - writer.add_timestep( 1, 5 * day, cfg.es, cfg.regionCache, cfg.wells, cfg.solution); - writer.add_timestep( 2, 10 * day, cfg.es, cfg.regionCache, cfg.wells, cfg.solution); + writer.add_timestep( 1, 2 * day, cfg.es, cfg.wells, cfg.solution); + writer.add_timestep( 1, 5 * day, cfg.es, cfg.wells, cfg.solution); + writer.add_timestep( 2, 10 * day, cfg.es, cfg.wells, cfg.solution); writer.write(); } @@ -729,9 +727,9 @@ BOOST_AUTO_TEST_CASE(region_production) { { out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); } @@ -757,9 +755,9 @@ BOOST_AUTO_TEST_CASE(region_injection) { setup cfg( "region_injection" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); @@ -786,9 +784,9 @@ BOOST_AUTO_TEST_CASE(BLOCK_VARIABLES) { setup cfg( "region_injection" ); out::Summary writer( cfg.es, cfg.config, cfg.grid, cfg.name ); - writer.add_timestep( 0, 0 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 1, 1 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); - writer.add_timestep( 2, 2 * day, cfg.es, cfg.regionCache, cfg.wells , cfg.solution); + writer.add_timestep( 0, 0 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 1, 1 * day, cfg.es, cfg.wells , cfg.solution); + writer.add_timestep( 2, 2 * day, cfg.es, cfg.wells , cfg.solution); writer.write(); auto res = readsum( cfg.name ); diff --git a/tests/test_writenumwells.cpp b/tests/test_writenumwells.cpp index c18ea4465..b215ea2c3 100644 --- a/tests/test_writenumwells.cpp +++ b/tests/test_writenumwells.cpp @@ -27,7 +27,7 @@ #include #include -#include +#include #include #include #include @@ -140,7 +140,7 @@ BOOST_AUTO_TEST_CASE(EclipseWriteRestartWellInfo) { auto es = Parser::parse( eclipse_data_filename, ParseContext() ); const auto num_cells = es.getInputGrid().getCartesianSize(); - EclipseWriter eclipseWriter( es, es.getInputGrid() ); + EclipseIO eclipseWriter( es, es.getInputGrid() ); int countTimeStep = es.getSchedule().getTimeMap().numTimesteps();