diff --git a/opm/output/Output.cpp b/opm/output/Output.cpp deleted file mode 100644 index e69de29bb..000000000 diff --git a/opm/output/Output.hpp b/opm/output/Output.hpp deleted file mode 100644 index e69de29bb..000000000 diff --git a/opm/output/OutputWriter.cpp b/opm/output/OutputWriter.cpp new file mode 100644 index 000000000..b8368e9ed --- /dev/null +++ b/opm/output/OutputWriter.cpp @@ -0,0 +1,103 @@ +#include "OutputWriter.hpp" + +#include +#include +#include +#include + +#include +#include +#include // unique_ptr + +using namespace std; +using namespace Opm; +using namespace Opm::parameter; + +namespace { + +/// Multiplexer over a list of output writers +struct MultiWriter : public OutputWriter { + /// Shorthand for a list of owned output writers + typedef forward_list > writers_t; + typedef writers_t::iterator it_t; + typedef unique_ptr ptr_t; + + /// Adopt a list of writers + MultiWriter (ptr_t writers) : writers_ (std::move (writers)) { } + + /// Forward the call to all writers + virtual void writeInit(const SimulatorTimerInterface &timer) { + for (it_t it = writers_->begin (); it != writers_->end (); ++it) { + (*it)->writeInit (timer); + } + } + + virtual void writeTimeStep(const SimulatorTimerInterface& timer, + const SimulatorState& reservoirState, + const WellState& wellState, + bool isSubstep) { + for (it_t it = writers_->begin (); it != writers_->end(); ++it) { + (*it)->writeTimeStep (timer, reservoirState, wellState, isSubstep); + } + } + +private: + ptr_t writers_; +}; + +/// Psuedo-constructor, can appear in template +template unique_ptr +create (const ParameterGroup& params, + std::shared_ptr eclipseState, + const Opm::PhaseUsage &phaseUsage, + std::shared_ptr grid) { + return unique_ptr (new Format (params, + eclipseState, + phaseUsage, + grid->number_of_cells, + grid->global_cell)); +} + +/// Map between keyword in configuration and the corresponding +/// constructor function (type) that should be called when detected. +/// The writer must have a constructor which takes params and parser. +/// +/// If you want to add more possible writer formats, just add them +/// to the list below! +typedef map (*)( + const ParameterGroup&, + std::shared_ptr eclipseState, + const Opm::PhaseUsage &phaseUsage, + std::shared_ptr )> map_t; +map_t FORMATS = { + { "output_ecl", &create }, +}; + +} // anonymous namespace + +unique_ptr +OutputWriter::create (const ParameterGroup& params, + std::shared_ptr eclipseState, + const Opm::PhaseUsage &phaseUsage, + std::shared_ptr grid) { + // allocate a list which will be filled with writers. this list + // is initially empty (no output). + MultiWriter::ptr_t list (new MultiWriter::writers_t ()); + + // loop through the map and see if we can find the key that is + // specified there + typedef map_t::iterator map_it_t; + for (map_it_t it = FORMATS.begin (); it != FORMATS.end(); ++it) { + // keyword which would indicate that this format should be used + const std::string name (it->first); + + // invoke the constructor for the type if we found the keyword + // and put the pointer to this writer onto the list + if (params.getDefault (name, false)) { + list->push_front (it->second (params, eclipseState, phaseUsage, grid)); + } + } + + // create a multiplexer from the list of formats we found + return unique_ptr (new MultiWriter (std::move (list))); +} diff --git a/opm/output/OutputWriter.hpp b/opm/output/OutputWriter.hpp new file mode 100644 index 000000000..8ae515772 --- /dev/null +++ b/opm/output/OutputWriter.hpp @@ -0,0 +1,118 @@ +/* + Copyright (c) 2013 Uni Research 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 OPM_OUTPUT_WRITER_HPP +#define OPM_OUTPUT_WRITER_HPP + +#include // unique_ptr, shared_ptr +#include + +struct UnstructuredGrid; + +namespace Opm { + +// forward declaration +class EclipseState; +namespace parameter { class ParameterGroup; } +class SimulatorState; +class WellState; +struct PhaseUsage; + +/*! + * Interface for writing non-compositional (blackoil, two-phase) simulation + * state to files. + * + * Use the create() function to setup a chain of writer based on the + * configuration values, e.g. + * + * \example + * \code{.cpp} + * ParameterGroup params (argc, argv, false); + * auto parser = std::make_shared ( + * params.get ("deck_filename")); + * + * std::unique_ptr writer = + * OutputWriter::create (params, parser); + * + * // before the first timestep + * writer->writeInit (timer); + * + * // after each timestep + * writer->writeTimeStep (timer, state, wellState); + * + * \endcode + */ +class OutputWriter { +public: + /// Allow derived classes to be used in the unique_ptr that is returned + /// from the create() method. (Every class that should be delete'd should + /// have a proper constructor, and if the base class isn't virtual then + /// the compiler won't call the right one when the unique_ptr goes out of + /// scope). + virtual ~OutputWriter () { } + + /** + * Write the static data (grid, PVT curves, etc) to disk. + * + * This routine should be called before the first timestep (i.e. when + * timer.currentStepNum () == 0) + */ + virtual void writeInit(const SimulatorTimerInterface &timer) = 0; + + /*! + * \brief Write a blackoil reservoir state to disk for later inspection with + * visualization tools like ResInsight + * + * \param[in] timer The timer providing time, time step, etc. information + * \param[in] reservoirState The thermodynamic state of the reservoir + * \param[in] wellState The production/injection data for all wells + * + * This routine should be called after the timestep has been advanced, + * i.e. timer.currentStepNum () > 0. + */ + virtual void writeTimeStep(const SimulatorTimerInterface& timer, + const SimulatorState& reservoirState, + const WellState& wellState, + bool isSubstep) = 0; + + /*! + * Create a suitable set of output formats based on configuration. + * + * @param params Configuration properties. This function will setup a + * multiplexer of applicable output formats based on the + * desired configuration values. + * + * @param deck Input deck used to set up the simulation. + * + * @param eclipseState The internalized input deck. + * + * @return Pointer to a multiplexer to all applicable output formats. + * + * @see Opm::share_obj + */ + static std::unique_ptr + create (const parameter::ParameterGroup& params, + std::shared_ptr eclipseState, + const Opm::PhaseUsage &phaseUsage, + std::shared_ptr grid); +}; + +} // namespace Opm + +#endif /* OPM_OUTPUT_WRITER_HPP */ diff --git a/opm/output/eclipse/CornerpointChopper.hpp b/opm/output/eclipse/CornerpointChopper.hpp new file mode 100644 index 000000000..8cf5e9f03 --- /dev/null +++ b/opm/output/eclipse/CornerpointChopper.hpp @@ -0,0 +1,457 @@ +/* + Copyright 2010 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 . +*/ + +#ifndef OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED +#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ + + class CornerPointChopper + { + public: + CornerPointChopper(const std::string& file) + { + Opm::ParseMode parseMode; + Opm::ParserPtr parser(new Opm::Parser()); + deck_ = parser->parseFile(file , parseMode); + + metricUnits_.reset(Opm::UnitSystem::newMETRIC()); + + const auto& specgridRecord = deck_->getKeyword("SPECGRID").getRecord(0); + dims_[0] = specgridRecord.getItem("NX").get< int >(0); + dims_[1] = specgridRecord.getItem("NY").get< int >(0); + dims_[2] = specgridRecord.getItem("NZ").get< int >(0); + + int layersz = 8*dims_[0]*dims_[1]; + const std::vector& ZCORN = deck_->getKeyword("ZCORN").getRawDoubleData(); + botmax_ = *std::max_element(ZCORN.begin(), ZCORN.begin() + layersz/2); + topmin_ = *std::min_element(ZCORN.begin() + dims_[2]*layersz - layersz/2, + ZCORN.begin() + dims_[2]*layersz); + + abszmax_ = *std::max_element(ZCORN.begin(), ZCORN.end()); + abszmin_ = *std::min_element(ZCORN.begin(), ZCORN.end()); + + std::cout << "Parsed grdecl file with dimensions (" + << dims_[0] << ", " << dims_[1] << ", " << dims_[2] << ")" << std::endl; + } + + + + + const int* dimensions() const + { + return dims_; + } + + + + + const int* newDimensions() const + { + return new_dims_; + } + + + + + const std::pair zLimits() const + { + return std::make_pair(botmax_, topmin_); + } + + const std::pair abszLimits() const + { + return std::make_pair(abszmin_, abszmax_); + } + + + void verifyInscribedShoebox(int imin, int ilen, int imax, + int jmin, int jlen, int jmax, + double zmin, double zlen, double zmax) + { + if (imin < 0) { + std::cerr << "Error! imin < 0 (imin = " << imin << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (ilen > dims_[0]) { + std::cerr << "Error! ilen larger than grid (ilen = " << ilen <<")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (imax > dims_[0]) { + std::cerr << "Error! imax larger than input grid (imax = " << imax << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (jmin < 0) { + std::cerr << "Error! jmin < 0 (jmin = " << jmin << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (jlen > dims_[1]) { + std::cerr << "Error! jlen larger than grid (jlen = " << jlen <<")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (jmax > dims_[1]) { + std::cerr << "Error! jmax larger than input grid (jmax = " << jmax << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (zmin < abszmin_) { + std::cerr << "Error! zmin ("<< zmin << ") less than minimum ZCORN value ("<< abszmin_ << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (zmax > abszmax_) { + std::cerr << "Error! zmax ("<< zmax << ") larger than maximal ZCORN value ("<< abszmax_ << ")\n"; + throw std::runtime_error("Inconsistent user input."); + } + if (zlen > (abszmax_ - abszmin_)) { + std::cerr << "Error! zlen ("<< zlen <<") larger than maximal ZCORN (" << abszmax_ << ") minus minimal ZCORN ("<< abszmin_ <<")\n"; + throw std::runtime_error("Inconsistent user input."); + } + } + + void chop(int imin, int imax, int jmin, int jmax, double zmin, double zmax, bool resettoorigin=true) + { + new_dims_[0] = imax - imin; + new_dims_[1] = jmax - jmin; + + // Filter the coord field + const std::vector& COORD = deck_->getKeyword("COORD").getRawDoubleData(); + int num_coord = COORD.size(); + if (num_coord != 6*(dims_[0] + 1)*(dims_[1] + 1)) { + std::cerr << "Error! COORD size (" << COORD.size() << ") not consistent with SPECGRID\n"; + throw std::runtime_error("Inconsistent COORD and SPECGRID."); + } + int num_new_coord = 6*(new_dims_[0] + 1)*(new_dims_[1] + 1); + double x_correction = COORD[6*((dims_[0] + 1)*jmin + imin)]; + double y_correction = COORD[6*((dims_[0] + 1)*jmin + imin) + 1]; + new_COORD_.resize(num_new_coord, 1e100); + for (int j = jmin; j < jmax + 1; ++j) { + for (int i = imin; i < imax + 1; ++i) { + int pos = (dims_[0] + 1)*j + i; + int new_pos = (new_dims_[0] + 1)*(j-jmin) + (i-imin); + // Copy all 6 coordinates for a pillar. + std::copy(COORD.begin() + 6*pos, COORD.begin() + 6*(pos + 1), new_COORD_.begin() + 6*new_pos); + if (resettoorigin) { + // Substract lowest x value from all X-coords, similarly for y, and truncate in z-direction + new_COORD_[6*new_pos] -= x_correction; + new_COORD_[6*new_pos + 1] -= y_correction; + new_COORD_[6*new_pos + 2] = 0; + new_COORD_[6*new_pos + 3] -= x_correction; + new_COORD_[6*new_pos + 4] -= y_correction; + new_COORD_[6*new_pos + 5] = zmax-zmin; + } + } + } + + // Get the z limits, check if they must be changed to make a shoe-box. + // This means that zmin must be greater than or equal to the highest + // coordinate of the bottom surface, while zmax must be less than or + // equal to the lowest coordinate of the top surface. + int layersz = 8*dims_[0]*dims_[1]; + const std::vector& ZCORN = deck_->getKeyword("ZCORN").getRawDoubleData(); + int num_zcorn = ZCORN.size(); + if (num_zcorn != layersz*dims_[2]) { + std::cerr << "Error! ZCORN size (" << ZCORN.size() << ") not consistent with SPECGRID\n"; + throw std::runtime_error("Inconsistent ZCORN and SPECGRID."); + } + + zmin = std::max(zmin, botmax_); + zmax = std::min(zmax, topmin_); + if (zmin >= zmax) { + std::cerr << "Error: zmin >= zmax (zmin = " << zmin << ", zmax = " << zmax << ")\n"; + throw std::runtime_error("zmin >= zmax"); + } + std::cout << "Chopping subsample, i: (" << imin << "--" << imax << ") j: (" << jmin << "--" << jmax << ") z: (" << zmin << "--" << zmax << ")" << std::endl; + + // We must find the maximum and minimum k value for the given z limits. + // First, find the first layer with a z-coordinate strictly above zmin. + int kmin = -1; + for (int k = 0; k < dims_[2]; ++k) { + double layer_max = *std::max_element(ZCORN.begin() + k*layersz, ZCORN.begin() + (k + 1)*layersz); + if (layer_max > zmin) { + kmin = k; + break; + } + } + // Then, find the last layer with a z-coordinate strictly below zmax. + int kmax = -1; + for (int k = dims_[2]; k > 0; --k) { + double layer_min = *std::min_element(ZCORN.begin() + (k - 1)*layersz, ZCORN.begin() + k*layersz); + if (layer_min < zmax) { + kmax = k; + break; + } + } + new_dims_[2] = kmax - kmin; + + // Filter the ZCORN field, build mapping from new to old cells. + double z_origin_correction = 0.0; + if (resettoorigin) { + z_origin_correction = zmin; + } + new_ZCORN_.resize(8*new_dims_[0]*new_dims_[1]*new_dims_[2], 1e100); + new_to_old_cell_.resize(new_dims_[0]*new_dims_[1]*new_dims_[2], -1); + int cellcount = 0; + int delta[3] = { 1, 2*dims_[0], 4*dims_[0]*dims_[1] }; + int new_delta[3] = { 1, 2*new_dims_[0], 4*new_dims_[0]*new_dims_[1] }; + for (int k = kmin; k < kmax; ++k) { + for (int j = jmin; j < jmax; ++j) { + for (int i = imin; i < imax; ++i) { + new_to_old_cell_[cellcount++] = dims_[0]*dims_[1]*k + dims_[0]*j + i; + int old_ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]); + int new_ix = 2*((i-imin)*new_delta[0] + (j-jmin)*new_delta[1] + (k-kmin)*new_delta[2]); + int old_indices[8] = { old_ix, old_ix + delta[0], + old_ix + delta[1], old_ix + delta[1] + delta[0], + old_ix + delta[2], old_ix + delta[2] + delta[0], + old_ix + delta[2] + delta[1], old_ix + delta[2] + delta[1] + delta[0] }; + int new_indices[8] = { new_ix, new_ix + new_delta[0], + new_ix + new_delta[1], new_ix + new_delta[1] + new_delta[0], + new_ix + new_delta[2], new_ix + new_delta[2] + new_delta[0], + new_ix + new_delta[2] + new_delta[1], new_ix + new_delta[2] + new_delta[1] + new_delta[0] }; + for (int cc = 0; cc < 8; ++cc) { + new_ZCORN_[new_indices[cc]] = std::min(zmax, std::max(zmin, ZCORN[old_indices[cc]])) - z_origin_correction; + } + } + } + } + + filterIntegerField("ACTNUM", new_ACTNUM_); + filterDoubleField("PORO", new_PORO_); + filterDoubleField("NTG", new_NTG_); + filterDoubleField("SWCR", new_SWCR_); + filterDoubleField("SOWCR", new_SOWCR_); + filterDoubleField("PERMX", new_PERMX_); + filterDoubleField("PERMY", new_PERMY_); + filterDoubleField("PERMZ", new_PERMZ_); + filterIntegerField("SATNUM", new_SATNUM_); + } + + /// Return a sub-deck with fields corresponding to the selected subset. + Opm::DeckConstPtr subDeck() + { + Opm::DeckPtr subDeck(new Opm::Deck); + + Opm::DeckKeyword specGridKw("SPECGRID"); + Opm::DeckRecord specGridRecord; + + auto nxItem = Opm::DeckItem::make< int >("NX"); + auto nyItem = Opm::DeckItem::make< int >("NY"); + auto nzItem = Opm::DeckItem::make< int >("NZ"); + auto numresItem = Opm::DeckItem::make< int >("NUMRES"); + auto coordTypeItem = Opm::DeckItem::make< std::string >("COORD_TYPE"); + + nxItem.push_back(new_dims_[0]); + nyItem.push_back(new_dims_[1]); + nzItem.push_back(new_dims_[2]); + numresItem.push_back(1); + coordTypeItem.push_back("F"); + + specGridRecord.addItem(std::move(nxItem)); + specGridRecord.addItem(std::move(nyItem)); + specGridRecord.addItem(std::move(nzItem)); + specGridRecord.addItem(std::move(numresItem)); + specGridRecord.addItem(std::move(coordTypeItem)); + + specGridKw.addRecord(std::move(specGridRecord)); + + subDeck->addKeyword(std::move(specGridKw)); + addDoubleKeyword_(subDeck, "COORD", /*dimension=*/"Length", new_COORD_); + addDoubleKeyword_(subDeck, "ZCORN", /*dimension=*/"Length", new_ZCORN_); + addIntKeyword_(subDeck, "ACTNUM", new_ACTNUM_); + addDoubleKeyword_(subDeck, "PORO", /*dimension=*/"1", new_PORO_); + addDoubleKeyword_(subDeck, "NTG", /*dimension=*/"1", new_NTG_); + addDoubleKeyword_(subDeck, "SWCR", /*dimension=*/"1", new_SWCR_); + addDoubleKeyword_(subDeck, "SOWCR", /*dimension=*/"1", new_SOWCR_); + addDoubleKeyword_(subDeck, "PERMX", /*dimension=*/"Permeability", new_PERMX_); + addDoubleKeyword_(subDeck, "PERMY", /*dimension=*/"Permeability", new_PERMY_); + addDoubleKeyword_(subDeck, "PERMZ", /*dimension=*/"Permeability", new_PERMZ_); + addIntKeyword_(subDeck, "SATNUM", new_SATNUM_); + return subDeck; + } + void writeGrdecl(const std::string& filename) + { + // Output new versions of SPECGRID, COORD, ZCORN, ACTNUM, PERMX, PORO, SATNUM. + std::ofstream out(filename.c_str()); + if (!out) { + std::cerr << "Could not open file " << filename << "\n"; + throw std::runtime_error("Could not open output file."); + } + out << "SPECGRID\n" << new_dims_[0] << ' ' << new_dims_[1] << ' ' << new_dims_[2] + << " 1 F\n/\n\n"; + + out.precision(15); + out.setf(std::ios::scientific); + + outputField(out, new_COORD_, "COORD", /* nl = */ 3); + outputField(out, new_ZCORN_, "ZCORN", /* nl = */ 4); + outputField(out, new_ACTNUM_, "ACTNUM"); + outputField(out, new_PORO_, "PORO", 4); + if (hasNTG()) {outputField(out, new_NTG_, "NTG", 4);} + if (hasSWCR()) {outputField(out, new_SWCR_, "SWCR", 4);} + if (hasSOWCR()) {outputField(out, new_SOWCR_, "SOWCR", 4);} + outputField(out, new_PERMX_, "PERMX", 4); + outputField(out, new_PERMY_, "PERMY", 4); + outputField(out, new_PERMZ_, "PERMZ", 4); + outputField(out, new_SATNUM_, "SATNUM"); + } + bool hasNTG() const {return !new_NTG_.empty(); } + bool hasSWCR() const {return !new_SWCR_.empty(); } + bool hasSOWCR() const {return !new_SOWCR_.empty(); } + + private: + Opm::DeckConstPtr deck_; + std::shared_ptr metricUnits_; + + double botmax_; + double topmin_; + double abszmin_; + double abszmax_; + std::vector new_COORD_; + std::vector new_ZCORN_; + std::vector new_ACTNUM_; + std::vector new_PORO_; + std::vector new_NTG_; + std::vector new_SWCR_; + std::vector new_SOWCR_; + std::vector new_PERMX_; + std::vector new_PERMY_; + std::vector new_PERMZ_; + std::vector new_SATNUM_; + int dims_[3]; + int new_dims_[3]; + std::vector new_to_old_cell_; + + void addDoubleKeyword_(Opm::DeckPtr subDeck, + const std::string& keywordName, + const std::string& dimensionString, + const std::vector& data) + { + if (data.empty()) + return; + + Opm::DeckKeyword dataKw(keywordName); + Opm::DeckRecord dataRecord; + auto dataItem = Opm::DeckItem::make< double >("DATA"); + + for (size_t i = 0; i < data.size(); ++i) { + dataItem.push_back(data[i]); + } + + std::shared_ptr dimension = metricUnits_->parse(dimensionString); + dataItem.push_backDimension(/*active=*/dimension, /*default=*/dimension); + + dataRecord.addItem(std::move(dataItem)); + dataKw.addRecord(std::move(dataRecord)); + subDeck->addKeyword(std::move(dataKw)); + } + + void addIntKeyword_(Opm::DeckPtr subDeck, + const std::string& keywordName, + const std::vector& data) + { + if (data.empty()) + return; + + Opm::DeckKeyword dataKw(keywordName); + Opm::DeckRecord dataRecord; + auto dataItem = Opm::DeckItem::make< int >("DATA"); + + for (size_t i = 0; i < data.size(); ++i) { + dataItem.push_back(data[i]); + } + + dataRecord.addItem(std::move(dataItem)); + dataKw.addRecord(std::move(dataRecord)); + subDeck->addKeyword(std::move(dataKw)); + } + + template + void outputField(std::ostream& os, + const std::vector& field, + const std::string& keyword, + const typename std::vector::size_type nl = 20) + { + if (field.empty()) return; + + os << keyword << '\n'; + + typedef typename std::vector::size_type sz_t; + + const sz_t n = field.size(); + for (sz_t i = 0; i < n; ++i) { + os << field[i] + << (((i + 1) % nl == 0) ? '\n' : ' '); + } + if (n % nl != 0) { + os << '\n'; + } + os << "/\n\n"; + } + + + + template + void filterField(const std::vector& field, + std::vector& output_field) + { + int sz = new_to_old_cell_.size(); + output_field.resize(sz); + for (int i = 0; i < sz; ++i) { + output_field[i] = field[new_to_old_cell_[i]]; + } + } + + void filterDoubleField(const std::string& keyword, std::vector& output_field) + { + if (deck_->hasKeyword(keyword)) { + const std::vector& field = deck_->getKeyword(keyword).getRawDoubleData(); + filterField(field, output_field); + } + } + + void filterIntegerField(const std::string& keyword, std::vector& output_field) + { + if (deck_->hasKeyword(keyword)) { + const std::vector& field = deck_->getKeyword(keyword).getIntData(); + filterField(field, output_field); + } + } + + }; + +} + + + + +#endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED diff --git a/opm/output/eclipse/EclipseGridInspector.hpp b/opm/output/eclipse/EclipseGridInspector.hpp new file mode 100644 index 000000000..df0dba0c9 --- /dev/null +++ b/opm/output/eclipse/EclipseGridInspector.hpp @@ -0,0 +1,103 @@ +//=========================================================================== +// +// File: EclipseGridInspector.h +// +// Created: Mon Jun 2 09:46:08 2008 +// +// Author: Atgeirr F Rasmussen +// +// $Date$ +// +// Revision: $Id: EclipseGridInspector.h,v 1.2 2008/08/18 14:16:12 atgeirr Exp $ +// +//=========================================================================== + +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. + + 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_ECLIPSEGRIDINSPECTOR_HEADER +#define OPM_ECLIPSEGRIDINSPECTOR_HEADER + +#include +#include + +#include + +namespace Opm +{ + +/** + @brief A class for inspecting the contents of an eclipse file. + + Given an Eclipse deck, this class may be used to answer certain + queries about its contents. + + @author Atgeirr F. Rasmussen + @date 2008/06/02 09:46:08 +*/ +class EclipseGridInspector +{ +public: + /// Constructor taking a parser as argument. + /// The parser must already have read an Eclipse file. + EclipseGridInspector(Opm::DeckConstPtr deck); + + /// Assuming that the pillars are vertical, compute the + /// volume of the cell given by logical coordinates (i, j, k). + double cellVolumeVerticalPillars(int i, int j, int k) const; + + /// Assuming that the pillars are vertical, compute the + /// volume of the cell given by the cell index + double cellVolumeVerticalPillars(int cell_idx) const; + + /// Compute the average dip in x- and y-direction of the + /// cell tops and bottoms relative to the xy-plane + std::pair cellDips(int i, int j, int k) const; + std::pair cellDips(int cell_idx) const; + + // Convert global cell index to logical ijk-coordinates + std::array cellIdxToLogicalCoords(int cell_idx) const; + + /// Returns a vector with the outer limits of grid (in the grid's unit). + /// The vector contains [xmin, xmax, ymin, ymax, zmin, zmax], as + /// read from COORDS and ZCORN + std::array getGridLimits() const; + + /// Returns the extent of the logical cartesian grid + /// as number of cells in the (i, j, k) directions. + std::array gridSize() const; + + /// Returns the eight z-values associated with a given cell. + /// The ordering is such that i runs fastest. That is, with + /// L = low and H = high: + /// {LLL, HLL, LHL, HHL, LLH, HLH, LHH, HHH }. + std::array cellZvals(int i, int j, int k) const; + +private: + Opm::DeckConstPtr deck_; + int logical_gridsize_[3]; + void init_(); + void checkLogicalCoords(int i, int j, int k) const; +}; + +} // namespace Opm + +#endif // OPM_ECLIPSEGRIDINSPECTOR_HEADER + diff --git a/opm/output/eclipse/EclipseIOUtil.hpp b/opm/output/eclipse/EclipseIOUtil.hpp new file mode 100644 index 000000000..1b8d3c202 --- /dev/null +++ b/opm/output/eclipse/EclipseIOUtil.hpp @@ -0,0 +1,34 @@ +#ifndef ECLIPSE_IO_UTIL_HPP +#define ECLIPSE_IO_UTIL_HPP + +#include +#include +#include + +namespace Opm +{ +namespace EclipseIOUtil +{ + + template + void addToStripedData(const std::vector& data, std::vector& result, size_t offset, size_t stride) { + int dataindex = 0; + for (size_t index = offset; index < result.size(); index += stride) { + result[index] = data[dataindex]; + ++dataindex; + } + } + + + template + void extractFromStripedData(const std::vector& data, std::vector& result, size_t offset, size_t stride) { + for (size_t index = offset; index < data.size(); index += stride) { + result.push_back(data[index]); + } + } + + +} //namespace EclipseIOUtil +} //namespace Opm + +#endif //ECLIPSE_IO_UTIL_HPP diff --git a/opm/output/eclipse/EclipseReader.hpp b/opm/output/eclipse/EclipseReader.hpp new file mode 100644 index 000000000..70c46f7d1 --- /dev/null +++ b/opm/output/eclipse/EclipseReader.hpp @@ -0,0 +1,35 @@ +#ifndef ECLIPSEREADER_HPP +#define ECLIPSEREADER_HPP + +#include +#include +#include +#include +#include + +namespace Opm +{ +/// +/// \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 +/// + + void init_from_restart_file(EclipseStateConstPtr eclipse_state, + int numcells, + const PhaseUsage& pu, + SimulatorState& simulator_state, + WellState& wellstate); + + +} + +#endif // ECLIPSEREADER_HPP diff --git a/opm/output/eclipse/EclipseUnits.hpp b/opm/output/eclipse/EclipseUnits.hpp new file mode 100644 index 000000000..866ba0cd3 --- /dev/null +++ b/opm/output/eclipse/EclipseUnits.hpp @@ -0,0 +1,63 @@ +/* + Copyright 2010 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 . +*/ + +#ifndef OPM_ECLIPSEUNITS_HEADER_INCLUDED +#define OPM_ECLIPSEUNITS_HEADER_INCLUDED + +namespace Opm +{ + struct EclipseUnits + { + double length; + double time; + double density; + double polymer_density; + double pressure; + double compressibility; + double viscosity; + double permeability; + double liqvol_s; + double liqvol_r; + double gasvol_s; + double gasvol_r; + double transmissibility; + + void setToOne() + { + length = 1.0; + time = 1.0; + density = 1.0; + polymer_density = 1.0; + pressure = 1.0; + compressibility = 1.0; + viscosity = 1.0; + permeability = 1.0; + liqvol_s = 1.0; + liqvol_r = 1.0; + gasvol_s = 1.0; + gasvol_r = 1.0; + transmissibility = 1.0; + } + }; + + +} // namespace Opm + + +#endif // OPM_ECLIPSEUNITS_HEADER_INCLUDED diff --git a/opm/output/eclipse/EclipseWriteRFTHandler.hpp b/opm/output/eclipse/EclipseWriteRFTHandler.hpp new file mode 100644 index 000000000..b0d8f9ed6 --- /dev/null +++ b/opm/output/eclipse/EclipseWriteRFTHandler.hpp @@ -0,0 +1,79 @@ +/* + Copyright 2015 Statoil ASA. + + 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_ECLIPSE_WRITE_RFT_HANDLER_HPP +#define OPM_ECLIPSE_WRITE_RFT_HANDLER_HPP + +#include +#include +#include + +#include + +#include +#include + + +namespace Opm { + class EclipseGrid; + class Well; + +namespace EclipseWriterDetails { + + + class EclipseWriteRFTHandler { + + public: + EclipseWriteRFTHandler(const int * compressedToCartesianCellIdx, size_t numCells, size_t cartesianSize); + + + void writeTimeStep(const std::string& filename, + const ert_ecl_unit_enum ecl_unit, + const SimulatorTimerInterface& simulatorTimer, + std::vector>& wells, + std::shared_ptr< const EclipseGrid > eclipseGrid, + std::vector& pressure, + std::vector& swat, + std::vector& sgas); + + + + private: + + ecl_rft_node_type * createEclRFTNode(std::shared_ptr< const Well > well, + const SimulatorTimerInterface& simulatorTimer, + std::shared_ptr< const EclipseGrid > eclipseGrid, + const std::vector& pressure, + const std::vector& swat, + const std::vector& sgas); + + void initGlobalToActiveIndex(const int * compressedToCartesianCellIdx, size_t numCells, size_t cartesianSize); + + std::vector globalToActiveIndex_; + + }; + + + + +}//namespace EclipseWriterDetails +}//namespace Opm + + +#endif //OPM_ECLIPSE_WRITE_RFT_HANDLER_HPP diff --git a/opm/output/eclipse/EclipseWriter.hpp b/opm/output/eclipse/EclipseWriter.hpp new file mode 100644 index 000000000..19d1e83f4 --- /dev/null +++ b/opm/output/eclipse/EclipseWriter.hpp @@ -0,0 +1,139 @@ +/* + Copyright (c) 2013 Andreas Lauser + Copyright (c) 2013 Uni Research AS + Copyright (c) 2014 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 OPM_ECLIPSE_WRITER_HPP +#define OPM_ECLIPSE_WRITER_HPP + +#include +#include +#include // WellType +#include + +#include +#include +#include + +#include + +#include +#include +#include +#include + +struct UnstructuredGrid; + +namespace Opm { + +// forward declarations +namespace EclipseWriterDetails { +class Summary; +} + +class SimulatorState; +class WellState; + +namespace parameter { class ParameterGroup; } + +/*! + * \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 EclipseWriter : public OutputWriter +{ +public: + /*! + * \brief Sets the common attributes required to write eclipse + * binary files using ERT. + */ + EclipseWriter(const parameter::ParameterGroup& params, + Opm::EclipseStateConstPtr eclipseState, + const Opm::PhaseUsage &phaseUsage, + int numCells, + const int* compressedToCartesianCellIdx); + + /** + * We need a destructor in the compilation unit to avoid the + * EclipseSummary being a complete type here. + */ + virtual ~EclipseWriter (); + + /** + * Write the static eclipse data (grid, PVT curves, etc) to disk. + */ + virtual void writeInit(const SimulatorTimerInterface &timer); + + /*! + * \brief Write a reservoir state and summary information to disk. + * + * + * The reservoir state can be inspected with visualization tools like + * ResInsight. + * + * The summary information can then be visualized using tools from + * ERT or ECLIPSE. Note that calling this method is only + * meaningful after the first time step has been completed. + * + * \param[in] timer The timer providing time step and time information + * \param[in] reservoirState The thermodynamic state of the reservoir + * \param[in] wellState The production/injection data for all wells + */ + virtual void writeTimeStep(const SimulatorTimerInterface& timer, + const SimulatorState& reservoirState, + const WellState& wellState, + bool isSubstep); + + + static int eclipseWellTypeMask(WellType wellType, WellInjector::TypeEnum injectorType); + static int eclipseWellStatusMask(WellCommon::StatusEnum wellStatus); + static ert_ecl_unit_enum convertUnitTypeErtEclUnitEnum(UnitSystem::UnitType unit); + +private: + Opm::EclipseStateConstPtr eclipseState_; + int numCells_; + std::array cartesianSize_; + const int* compressedToCartesianCellIdx_; + std::vector< int > gridToEclipseIdx_; + double deckToSiPressure_; + double deckToSiTemperatureFactor_; + double deckToSiTemperatureOffset_; + bool enableOutput_; + int writeStepIdx_; + int reportStepIdx_; + std::string outputDir_; + std::string baseName_; + PhaseUsage phaseUsage_; // active phases in the input deck + std::shared_ptr summary_; + + void init(const parameter::ParameterGroup& params); +}; + +typedef std::shared_ptr EclipseWriterPtr; +typedef std::shared_ptr EclipseWriterConstPtr; + +} // namespace Opm + + +#endif // OPM_ECLIPSE_WRITER_HPP diff --git a/opm/output/eclipse/writeECLData.hpp b/opm/output/eclipse/writeECLData.hpp new file mode 100644 index 000000000..e962573af --- /dev/null +++ b/opm/output/eclipse/writeECLData.hpp @@ -0,0 +1,46 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + Copyright 2012 Statoil ASA. + + 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_WRITEECLDATA_HEADER_INCLUDED +#define OPM_WRITEECLDATA_HEADER_INCLUDED + + +#include +#include + +#include + +struct UnstructuredGrid; + +namespace Opm +{ + + // ECLIPSE output for general grids. + void writeECLData(const UnstructuredGrid& grid, + const DataMap& data, + const int current_step, + const double current_time, + const boost::posix_time::ptime& current_date_time, + const std::string& output_dir, + const std::string& base_name); + +} + +#endif diff --git a/opm/output/vag/vag.cpp b/opm/output/vag/vag.cpp new file mode 100644 index 000000000..2e91ce4b6 --- /dev/null +++ b/opm/output/vag/vag.cpp @@ -0,0 +1,408 @@ +/*=========================================================================== +// +// File: vag.cpp +// +// Created: 2012-06-08 15:45:53+0200 +// +// Authors: Knut-Andreas Lie +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Xavier Raynaud +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + Copyright 2012 Statoil ASA. + + 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +namespace Opm +{ + void readPosStruct(std::istream& is,int n,PosStruct& pos_struct){ + using namespace std; + //PosStruct pos_struct; + pos_struct.pos.resize(n+1); + pos_struct.pos[0]=0; + for(int i=0;i< n;++i){ + int number; + is >> number ; + //cout <> value; + // cout << value << " "; + pos_struct.value.push_back(value); + } + //cout << endl; + } + if(int(pos_struct.value.size()) != pos_struct.pos[n]){ + cerr << "Failed to read pos structure" << endl; + cerr << "pos_struct.value.size()" << pos_struct.value.size() << endl; + cerr << "pos_struct.pos[n+1]" << pos_struct.pos[n] << endl; + } + } + void writePosStruct(std::ostream& os,PosStruct& pos_struct){ + using namespace std; + //PosStruct pos_struct; + if(pos_struct.pos.size()==0){ + return; + } + int n=pos_struct.pos.size()-1; + pos_struct.pos.resize(n+1); + pos_struct.pos[0]=0; + for(int i=0;i< n;++i){ + int number=pos_struct.pos[i+1]-pos_struct.pos[i]; + os << number << " "; + for(int j=0;j< number;++j){ + os << pos_struct.value[pos_struct.pos[i]+j] << " "; + } + os << endl; + } + } + void readVagGrid(std::istream& is,Opm::VAG& vag_grid){ + using namespace std; + using namespace Opm; + while (!is.eof()) { + string keyword; + is >> keyword; + //cout << keyword<< endl; + if(keyword == "Number"){ + string stmp; + is >> stmp; + if(stmp == "of"){ + string entity; + is >> entity; + getline(is,stmp); + int number; + is >> number; + if(entity=="vertices"){ + vag_grid.number_of_vertices=number; + }else if((entity=="volumes") || (entity=="control")){ + vag_grid.number_of_volumes=number; + }else if(entity=="faces"){ + vag_grid.number_of_faces=number; + }else if(entity=="edges"){ + vag_grid.number_of_edges=number; + } + cout << "Found Number of: " << entity <<" " << number << endl; + } else { + cerr << "Wrong format: Not of after Number" << endl; + return; + } + }else{ + // read geometry defined by vertices + if(keyword=="Vertices"){ + int number; + is >> number; + vag_grid.vertices.resize(3*number);// assume 3d data + readVector(is,vag_grid.vertices); + } + // here starts the reding of all pos structures + else if(keyword=="Volumes->Faces" || keyword=="Volumes->faces"){ + //vag_grid.volumes_to_faces= + int number; + is >> number; + readPosStruct(is,number,vag_grid.volumes_to_faces); + cout << "Volumes->Faces: Number of " << number << endl; + }else if(keyword=="Faces->edges" || keyword=="Faces->Edges" || keyword=="Faces->Edgess"){ + int number; + is >> number; + //vag_grid.volumes_to_faces= + readPosStruct(is,number,vag_grid.faces_to_edges); + cout << "Faces->edges: Number of " << number << endl; + }else if(keyword=="Faces->Vertices" || keyword=="Faces->vertices"){ + int number; + is >> number; + //vag_grid.volumes_to_faces= + readPosStruct(is,number,vag_grid.faces_to_vertices); + cout << "Faces->Vertices: Number of " << number << endl; + }else if(keyword=="Volumes->Vertices" || keyword=="Volumes->Verticess"){ + int number; + is >> number; + //vag_grid.volumes_to_faces= + readPosStruct(is,number,vag_grid.volumes_to_vertices); + cout << "Volumes->Vertices: Number of " << number << endl; + } + + // read simple mappings + else if(keyword=="Edge" || keyword=="Edges"){ + int number; + is >> number; + vag_grid.edges.resize(2*number); + readVector(is,vag_grid.edges); + cout << "Edges: Number of " << number << endl; + }else if(keyword=="Faces->Volumes" || keyword=="Faces->Control"){ + int number; + if(keyword=="Faces->Control"){ + string vol; + is >> vol; + } + is >> number; + vag_grid.faces_to_volumes.resize(2*number); + readVector(is,vag_grid.faces_to_volumes); + cout << "Faces->Volumes: Number of " << number << endl; + } + // read material + else if(keyword=="Material"){ + string snum; + is >> snum; + int number; + is >> number; + cout << "Material number " << number << endl; + // we read all the rest into doubles + while(!is.eof()){ + double value; + is >> value; + //cout << value << endl; + vag_grid.material.push_back(value); + } + }else{ + //cout << "keyword; + } + //cout << "Found" << keyword << "Number of " << number << endl; + } + } + } + void vagToUnstructuredGrid(Opm::VAG& vag_grid,UnstructuredGrid& grid){ + using namespace std; + using namespace Opm; + cout << "Converting grid" << endl; + cout << "Warning:: orignial grid may not be edge confomal" << endl; + cout << " inverse mappings from edges will be wrong" << endl; + grid.dimensions=3; + grid.number_of_cells=vag_grid.number_of_volumes; + grid.number_of_faces=vag_grid.number_of_faces; + grid.number_of_nodes=vag_grid.number_of_vertices; + + // fill face_nodes + for(int i=0;i< int(vag_grid.faces_to_vertices.pos.size());++i){ + grid.face_nodepos[i] = vag_grid.faces_to_vertices.pos[i]; + } + for(int i=0;i< int(vag_grid.faces_to_vertices.value.size());++i){ + grid.face_nodes[i] = vag_grid.faces_to_vertices.value[i]-1; + } + // fill cell_face + for(int i=0;i< int(vag_grid.volumes_to_faces.pos.size());++i){ + grid.cell_facepos[i] = vag_grid.volumes_to_faces.pos[i]; + } + for(int i=0;i< int(vag_grid.volumes_to_faces.value.size());++i){ + grid.cell_faces[i] = vag_grid.volumes_to_faces.value[i]-1; + } + // fill face_cells + for(int i=0;i< int(vag_grid.faces_to_volumes.size());++i){ + grid.face_cells[i] = vag_grid.faces_to_volumes[i]-1; + } + + // fill node_cordinates. This is the only geometry given in the vag + for(int i=0;i< int(vag_grid.vertices.size());++i){ + grid.node_coordinates[i] = vag_grid.vertices[i]; + } + // informations in edges, faces_to_eges, faces_to_vertices, volume_to_vertices and materials + // is not used + cout << "Computing geometry" << endl; + compute_geometry(&grid); + + } + + void unstructuredGridToVag(UnstructuredGrid& grid,Opm::VAG& vag_grid){ + using namespace std; + using namespace Opm; + cout << "Converting grid" << endl; + // grid.dimensions=3; + vag_grid.number_of_volumes=grid.number_of_cells; + vag_grid.number_of_faces=grid.number_of_faces; + vag_grid.number_of_vertices=grid.number_of_nodes; + + // resizing vectors + vag_grid.vertices.resize(grid.number_of_nodes*3); + vag_grid.faces_to_vertices.pos.resize(grid.number_of_faces+1); + vag_grid.faces_to_vertices.value.resize(grid.face_nodepos[grid.number_of_faces]); + vag_grid.faces_to_volumes.resize(2*grid.number_of_faces); + vag_grid.volumes_to_faces.pos.resize(grid.number_of_cells+1); + vag_grid.volumes_to_faces.value.resize(grid.cell_facepos[grid.number_of_cells]);//not known + + + + + // fill face_nodes + for(int i=0;i< int(vag_grid.faces_to_vertices.pos.size());++i){ + vag_grid.faces_to_vertices.pos[i] = grid.face_nodepos[i]; + } + + for(int i=0;i< int(vag_grid.faces_to_vertices.value.size());++i){ + vag_grid.faces_to_vertices.value[i] = grid.face_nodes[i] +1; + } + + // fill cell_face + for(int i=0;i< int(vag_grid.volumes_to_faces.pos.size());++i){ + vag_grid.volumes_to_faces.pos[i] = grid.cell_facepos[i]; + } + for(int i=0;i< int(vag_grid.volumes_to_faces.value.size());++i){ + vag_grid.volumes_to_faces.value[i] = grid.cell_faces[i] +1; + } + // fill face_cells + for(int i=0;i< int(vag_grid.faces_to_volumes.size());++i){ + vag_grid.faces_to_volumes[i] = grid.face_cells[i] +1; + } + + // fill node_cordinates. This is the only geometry given in the vag + for(int i=0;i< int(vag_grid.vertices.size());++i){ + vag_grid.vertices[i] = grid.node_coordinates[i]; + } + + + // The missing field need to be constructed + // gennerate volume to vertice mapping + std::vector< std::set > volumes_to_vertices(grid.number_of_cells); + for(int i=0;i < grid.number_of_cells; ++i){ + int nlf=grid.cell_facepos[i+1]-grid.cell_facepos[i]; + std::set nodes; + for(int j=0; j < nlf; ++j){ + int face = grid.cell_faces[grid.cell_facepos[i]+j]; + int nlv = grid.face_nodepos[face+1]-grid.face_nodepos[face]; + for(int k=0; k< nlv; ++k){ + int node = grid.face_nodes[grid.face_nodepos[face]+k]+1; + nodes.insert(node); + } + } + volumes_to_vertices[i]=nodes; + } + + // fill volume to vertice map + vag_grid.volumes_to_vertices.pos.resize(grid.number_of_cells+1); + vag_grid.volumes_to_vertices.value.resize(0); + vag_grid.volumes_to_vertices.pos[0]=0; + for(int i=0;i < grid.number_of_cells;++i){ + int nv=volumes_to_vertices[i].size(); + vag_grid.volumes_to_vertices.pos[i+1]=vag_grid.volumes_to_vertices.pos[i]+nv; + std::set::iterator it; + for(it=volumes_to_vertices[i].begin();it!=volumes_to_vertices[i].end();++it){ + vag_grid.volumes_to_vertices.value.push_back(*it); + } + } + + std::set< std::set > edges; + std::vector< std::vector< std::set > > faces_spares; + int nfe=0; + faces_spares.resize(grid.number_of_faces); + for(int i=0;i < grid.number_of_faces;++i){ + int ne=grid.face_nodepos[i+1]-grid.face_nodepos[i]; + nfe=nfe+ne; + + for(int j=0; j < ne-1;++j){ + int node1=grid.face_nodes[grid.face_nodepos[i]+j]+1; + int node2=grid.face_nodes[grid.face_nodepos[i]+j+1]+1; + std::set spair; + spair.insert(node1); + spair.insert(node2); + edges.insert(spair); + faces_spares[i].push_back(spair); + } + // add end segment + { + std::set spair; + int node1=grid.face_nodes[grid.face_nodepos[i]+ne-1]+1; + int node2=grid.face_nodes[grid.face_nodepos[i]]+1; + spair.insert(node1); + spair.insert(node2); + edges.insert(spair); + faces_spares[i].push_back(spair); + } + } + + // make edge numbering and fill edges + std::map, int> edge_map; + std::set< std::set >::iterator it; + vag_grid.edges.resize(0); + int k=0; + for(it=edges.begin(); it!=edges.end();++it){ + edge_map.insert(std::pair< std::set , int >(*it,k)); + k=k+1; + std::set::iterator sit; + for(sit=(*it).begin();sit!=(*it).end();++sit){ + vag_grid.edges.push_back(*sit); + } + } + // fill face_to_egdes + vag_grid.number_of_edges=edges.size(); + vag_grid.faces_to_edges.pos.resize(vag_grid.number_of_faces+1); + for(int i=0;i < grid.number_of_faces;++i){ + int ne=grid.face_nodepos[i+1]-grid.face_nodepos[i]; + vag_grid.faces_to_edges.pos[i+1]=vag_grid.faces_to_edges.pos[i]+ne; + for(int j=0;jfaces " << vag_grid.volumes_to_faces.pos.size()-1 << endl; + writePosStruct(os, vag_grid.volumes_to_faces); + os << "Volumes->Vertices " << vag_grid.volumes_to_vertices.pos.size()-1 << endl; + writePosStruct(os, vag_grid.volumes_to_vertices); + os << "Faces->edges " << vag_grid.faces_to_edges.pos.size()-1 << endl; + writePosStruct(os, vag_grid.faces_to_edges); + os << "Faces->vertices " << vag_grid.faces_to_vertices.pos.size()-1 << endl; + writePosStruct(os, vag_grid.faces_to_vertices); + os << "Faces->Control volumes " << floor(vag_grid.faces_to_volumes.size()/2) << endl; + writeVector(os,vag_grid.faces_to_volumes,2); + os << "Edges " << floor(vag_grid.edges.size()/2) << endl; + writeVector(os,vag_grid.edges,2); + /* + assert(vag_grid.material.size()%vag_grid.number_of_volumes==0); + int lines= floor(vag_grid.material.size()/vag_grid.number_of_volumes); + os << "Material number " << 1 << endl; + writeVector(os,vag_grid.material,lines); + */ + + } + + +} diff --git a/opm/output/vag/vag.hpp b/opm/output/vag/vag.hpp new file mode 100644 index 000000000..ab84af774 --- /dev/null +++ b/opm/output/vag/vag.hpp @@ -0,0 +1,165 @@ +/*=========================================================================== +// +// File: vag.hpp +// +// Created: 2012-06-08 15:46:23+0200 +// +// Authors: Knut-Andreas Lie +// Halvor M. Nilsen +// Atgeirr F. Rasmussen +// Xavier Raynaud +// Bård Skaflestad +// +//==========================================================================*/ + + +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + Copyright 2012 Statoil ASA. + + 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_VAG_HPP_HEADER +#define OPM_VAG_HPP_HEADER + + + +#include +#include +#include +#include +#include + +namespace Opm +{ + /** + Struct to hold mapping from the natural numbers less than pos.size()-1 to + a set of integers. value(pos(i):pos(i+1)-1) hold the integers corresponding to i. + pos(end)-1==value.size(); + */ + struct PosStruct{ + std::vector pos; + std::vector value; + }; + /** + Structure to represent the unstructured vag grid format. The format is only for + 3D grids. + */ + struct VAG{ + int number_of_vertices; + int number_of_volumes; + int number_of_faces; + int number_of_edges; + /** Vertices. The coordinates of vertice i is [vetices[3*i:3*i+2]*/ + std::vector vertices; + /** Mapping from volumes to faces */ + PosStruct volumes_to_faces; + /** Mapping from volumes to vertices */ + PosStruct volumes_to_vertices; + /** Mapping from faces to edges */ + PosStruct faces_to_edges; + /** Mapping from faces to vertices */ + PosStruct faces_to_vertices; + /** The edge i is given by the nodes edges[2*i:2*i+1] */ + std::vector edges; + /** The two neigbours of the face i is faces_to_volumes[2*i:2*i+1] */ + std::vector faces_to_volumes; + /** A vector containing information of each volume. The size is n*number_of_volumes. + For each i this is the information: + material[n*i] is the volume number and should be transformed to integer + material[n*i+1] is a tag and should be transformed to integer + material[n*i+2:n*(i+1)-1] represent propertices. + */ + std::vector material; + }; + /** + Function the vag grid format and make a vag_grid struct. This structure + is intended to be converted to a grid. + \param[in] is is is stream of the file. + \param[out] vag_grid is a reference to a vag_grid struct. + */ + void readVagGrid(std::istream& is,Opm::VAG& vag_grid); + /** + Function to write vag format. + \param[out] is is is stream of the file. + \param[in] vag_grid is a reference to a vag_grid struct. + + */ + void writeVagFormat(std::ostream& os,Opm::VAG& vag_grid); + /** + Function to read a vector of some type from a stream. + \param[in] os is is stream of the file. + \param[out] vag_grid is a resized and filled vector containing the quantiy read. + */ + template + void readVector(std::istream& is,std::vector& vec){ + using namespace std; + for(int i=0;i< int(vec.size());++i){ + is >> vec[i]; + } + } + /** + Function to write a vector of some type from a stream. + \param[in] os is is stream of the file. + \param[out] vag_grid is a resized and filled vector containing the quantiy read. + \param[in] n number of doubles on each line. + */ + template + void writeVector(std::ostream& os,std::vector& vec,int n){ + typedef typename std::vector::size_type sz_t; + + const sz_t nn = n; + + for (sz_t i = 0; i < vec.size(); ++i) { + os << vec[i] << (((i % nn) == 0) ? '\n' : ' '); + } + + if ((vec.size() % nn) != 0) { + os << '\n'; + } + } + + /** + Read pos struct type mapping from a stream + \param[in] is is stream + \param[in] n number of lines to read + \param[out] pos_struct reference to PosStruct + */ + void readPosStruct(std::istream& is,int n,PosStruct& pos_struct); + /** + Read pos struct type mapping from a stream + \param[in] os is stream to write to + \param[in] pos_struct to write + */ + void writePosStruct(std::ostream& os,PosStruct& pos_struct); + + /** + Fill a UnstructuredGrid from a vag_grid. + \param[out] vag_grid s is a valid vag_grid struct. + \param[in] grid is a grid with have allocated correct size to each pointer. + */ + void vagToUnstructuredGrid(Opm::VAG& vag_grid,UnstructuredGrid& grid); + + /** + Fill a vag_grid from UnstructuredGrid + \param[out] vag_grid s is a valid vag_grid struct. + \param[in] grid is a grid with have allocated correct size to each pointer. + */ + void unstructuredGridToVag(UnstructuredGrid& grid, Opm::VAG& vag_grid); +} +#endif /* OPM_VAG_HPP_HEADER */ + diff --git a/opm/output/vtk/writeVtkData.cpp b/opm/output/vtk/writeVtkData.cpp new file mode 100644 index 000000000..5d88bc7c8 --- /dev/null +++ b/opm/output/vtk/writeVtkData.cpp @@ -0,0 +1,319 @@ +/* + Copyright 2012 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 +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + + + +namespace Opm +{ + + void writeVtkData(const std::array& dims, + const std::array& cell_size, + const DataMap& data, + std::ostream& os) + { + // Dimension is hardcoded in the prototype and the next two lines, + // but the rest is flexible (allows dimension == 2 or 3). + int dimension = 3; + int num_cells = dims[0]*dims[1]*dims[2]; + + assert(dimension == 2 || dimension == 3); + assert(num_cells == dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2])); + + os << "# vtk DataFile Version 2.0\n"; + os << "Structured Grid\n \n"; + os << "ASCII \n"; + os << "DATASET STRUCTURED_POINTS\n"; + + os << "DIMENSIONS " + << dims[0] + 1 << " " + << dims[1] + 1 << " "; + if (dimension == 3) { + os << dims[2] + 1; + } else { + os << 1; + } + os << "\n"; + + os << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n"; + + os << "SPACING " << cell_size[0] << " " << cell_size[1]; + if (dimension == 3) { + os << " " << cell_size[2]; + } else { + os << " " << 0.0; + } + os << "\n"; + + os << "\nCELL_DATA " << num_cells << '\n'; + for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { + std::string name = dit->first; + os << "SCALARS " << name << " float" << '\n'; + os << "LOOKUP_TABLE " << name << "_table " << '\n'; + const std::vector& field = *(dit->second); + // We always print only the first data item for every + // cell, using 'stride'. + // This is a hack to get water saturation nicely. + // \TODO: Extend to properly printing vector data. + const int stride = field.size()/num_cells; + const int num_per_line = 5; + for (int c = 0; c < num_cells; ++c) { + os << field[stride*c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + } + + typedef std::map PMap; + + + struct Tag + { + Tag(const std::string& tag, const PMap& props, std::ostream& os) + : name_(tag), os_(os) + { + indent(os); + os << "<" << tag; + for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) { + os << " " << it->first << "=\"" << it->second << "\""; + } + os << ">\n"; + ++indent_; + } + Tag(const std::string& tag, std::ostream& os) + : name_(tag), os_(os) + { + indent(os); + os << "<" << tag << ">\n"; + ++indent_; + } + ~Tag() + { + --indent_; + indent(os_); + os_ << "\n"; + } + static void indent(std::ostream& os) + { + for (int i = 0; i < indent_; ++i) { + os << " "; + } + } + private: + static int indent_; + std::string name_; + std::ostream& os_; + }; + + int Tag::indent_ = 0; + + + void writeVtkData(const UnstructuredGrid& grid, + const DataMap& data, + std::ostream& os) + { + if (grid.dimensions != 3) { + OPM_THROW(std::runtime_error, "Vtk output for 3d grids only"); + } + os.precision(12); + os << "\n"; + PMap pm; + pm["type"] = "UnstructuredGrid"; + Tag vtkfiletag("VTKFile", pm, os); + Tag ugtag("UnstructuredGrid", os); + int num_pts = grid.number_of_nodes; + int num_cells = grid.number_of_cells; + pm.clear(); + pm["NumberOfPoints"] = boost::lexical_cast(num_pts); + pm["NumberOfCells"] = boost::lexical_cast(num_cells); + Tag piecetag("Piece", pm, os); + { + Tag pointstag("Points", os); + pm.clear(); + pm["type"] = "Float64"; + pm["Name"] = "Coordinates"; + pm["NumberOfComponents"] = "3"; + pm["format"] = "ascii"; + Tag datag("DataArray", pm, os); + for (int i = 0; i < num_pts; ++i) { + Tag::indent(os); + os << grid.node_coordinates[3*i + 0] << ' ' + << grid.node_coordinates[3*i + 1] << ' ' + << grid.node_coordinates[3*i + 2] << '\n'; + } + } + { + Tag cellstag("Cells", os); + pm.clear(); + pm["type"] = "Int32"; + pm["NumberOfComponents"] = "1"; + pm["format"] = "ascii"; + std::vector cell_numpts; + cell_numpts.reserve(num_cells); + { + pm["Name"] = "connectivity"; + Tag t("DataArray", pm, os); + int hf = 0; + for (int c = 0; c < num_cells; ++c) { + std::set cell_pts; + for (; hf < grid.cell_facepos[c+1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* fnbeg = grid.face_nodes + grid.face_nodepos[f]; + const int* fnend = grid.face_nodes + grid.face_nodepos[f+1]; + cell_pts.insert(fnbeg, fnend); + } + cell_numpts.push_back(cell_pts.size()); + Tag::indent(os); + std::copy(cell_pts.begin(), cell_pts.end(), + std::ostream_iterator(os, " ")); + os << '\n'; + } + } + { + pm["Name"] = "offsets"; + Tag t("DataArray", pm, os); + int offset = 0; + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + offset += cell_numpts[c]; + os << offset << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + std::vector cell_foffsets; + cell_foffsets.reserve(num_cells); + { + pm["Name"] = "faces"; + Tag t("DataArray", pm, os); + const int* fp = grid.cell_facepos; + int offset = 0; + for (int c = 0; c < num_cells; ++c) { + Tag::indent(os); + os << fp[c+1] - fp[c] << '\n'; + ++offset; + for (int hf = fp[c]; hf < fp[c+1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* np = grid.face_nodepos; + int f_num_pts = np[f+1] - np[f]; + Tag::indent(os); + os << f_num_pts << ' '; + ++offset; + std::copy(grid.face_nodes + np[f], + grid.face_nodes + np[f+1], + std::ostream_iterator(os, " ")); + os << '\n'; + offset += f_num_pts; + } + cell_foffsets.push_back(offset); + } + } + { + pm["Name"] = "faceoffsets"; + Tag t("DataArray", pm, os); + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << cell_foffsets[c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + { + pm["type"] = "UInt8"; + pm["Name"] = "types"; + Tag t("DataArray", pm, os); + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << "42 "; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + } + { + pm.clear(); + if (data.find("saturation") != data.end()) { + pm["Scalars"] = "saturation"; + } else if (data.find("pressure") != data.end()) { + pm["Scalars"] = "pressure"; + } + Tag celldatatag("CellData", pm, os); + pm.clear(); + pm["NumberOfComponents"] = "1"; + pm["format"] = "ascii"; + pm["type"] = "Float64"; + for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { + pm["Name"] = dit->first; + const std::vector& field = *(dit->second); + const int num_comps = field.size()/grid.number_of_cells; + pm["NumberOfComponents"] = boost::lexical_cast(num_comps); + Tag ptag("DataArray", pm, os); + const int num_per_line = num_comps == 1 ? 5 : num_comps; + for (int item = 0; item < num_cells*num_comps; ++item) { + if (item % num_per_line == 0) { + Tag::indent(os); + } + double value = field[item]; + if (std::fabs(value) < std::numeric_limits::min()) { + // Avoiding denormal numbers to work around + // bug in Paraview. + value = 0.0; + } + os << value << ' '; + if (item % num_per_line == num_per_line - 1 + || item == num_cells - 1) { + os << '\n'; + } + } + } + } + } + +} // namespace Opm + diff --git a/opm/output/vtk/writeVtkData.hpp b/opm/output/vtk/writeVtkData.hpp new file mode 100644 index 000000000..226750dd3 --- /dev/null +++ b/opm/output/vtk/writeVtkData.hpp @@ -0,0 +1,48 @@ +/* + Copyright 2012 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 . +*/ + +#ifndef OPM_WRITEVTKDATA_HEADER_INCLUDED +#define OPM_WRITEVTKDATA_HEADER_INCLUDED + + +#include +#include +#include +#include +#include +#include + +struct UnstructuredGrid; + +namespace Opm +{ + + /// Vtk output for cartesian grids. + void writeVtkData(const std::array& dims, + const std::array& cell_size, + const DataMap& data, + std::ostream& os); + + /// Vtk output for general grids. + void writeVtkData(const UnstructuredGrid& grid, + const DataMap& data, + std::ostream& os); +} // namespace Opm + +#endif // OPM_WRITEVTKDATA_HEADER_INCLUDED diff --git a/src/opm/output/eclipse/EclipseGridInspector.cpp b/src/opm/output/eclipse/EclipseGridInspector.cpp new file mode 100644 index 000000000..12d88498b --- /dev/null +++ b/src/opm/output/eclipse/EclipseGridInspector.cpp @@ -0,0 +1,350 @@ +//=========================================================================== +// +// File: EclipseGridInspector.C +// +// Created: Mon Jun 2 12:17:51 2008 +// +// Author: Atgeirr F Rasmussen +// +// $Date$ +// +// $Revision$ +// +// Revision: $Id: EclipseGridInspector.C,v 1.2 2008/08/18 14:16:13 atgeirr Exp $ +// +//=========================================================================== + +/* + Copyright 2009, 2010 SINTEF ICT, Applied Mathematics. + Copyright 2009, 2010 Statoil ASA. + + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +namespace Opm +{ +EclipseGridInspector::EclipseGridInspector(Opm::DeckConstPtr deck) + : deck_(deck) +{ + init_(); +} + +void EclipseGridInspector::init_() +{ + if (!deck_->hasKeyword("COORD")) { + OPM_THROW(std::runtime_error, "Needed field \"COORD\" is missing in file"); + } + if (!deck_->hasKeyword("ZCORN")) { + OPM_THROW(std::runtime_error, "Needed field \"ZCORN\" is missing in file"); + } + + if (deck_->hasKeyword("SPECGRID")) { + const auto& specgridRecord = + deck_->getKeyword("SPECGRID").getRecord(0); + logical_gridsize_[0] = specgridRecord.getItem("NX").get< int >(0); + logical_gridsize_[1] = specgridRecord.getItem("NY").get< int >(0); + logical_gridsize_[2] = specgridRecord.getItem("NZ").get< int >(0); + } else if (deck_->hasKeyword("DIMENS")) { + const auto& dimensRecord = + deck_->getKeyword("DIMENS").getRecord(0); + logical_gridsize_[0] = dimensRecord.getItem("NX").get< int >(0); + logical_gridsize_[1] = dimensRecord.getItem("NY").get< int >(0); + logical_gridsize_[2] = dimensRecord.getItem("NZ").get< int >(0); + } else { + OPM_THROW(std::runtime_error, "Found neither SPECGRID nor DIMENS in file. At least one is needed."); + } + +} + +/** + Return the dip slopes for the cell relative to xy-plane in x- and y- direction. + Dip slope is average rise in positive x-direction over cell length in x-direction. + Similarly for y. + + Current implementation is for vertical pillars, but is not difficult to fix. + + @returns a std::pair with x-dip in first component and y-dip in second. +*/ +std::pair EclipseGridInspector::cellDips(int i, int j, int k) const +{ + checkLogicalCoords(i, j, k); + const std::vector& pillc = + deck_->getKeyword("COORD").getSIDoubleData(); + int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1); + if (6*num_pillars != int(pillc.size())) { + throw std::runtime_error("Wrong size of COORD field."); + } + const std::vector& z = + deck_->getKeyword("ZCORN").getSIDoubleData(); + int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2]; + if (8*num_cells != int(z.size())) { + throw std::runtime_error("Wrong size of ZCORN field"); + } + + // Pick ZCORN-value for all 8 corners of the given cell + std::array cellz = cellZvals(i, j, k); + + // Compute rise in positive x-direction for all four edges (and then find mean) + // Current implementation is for regularly placed and vertical pillars! + int numxpill = logical_gridsize_[0] + 1; + int pix = i + j*numxpill; + double cell_xlength = pillc[6*(pix + 1)] - pillc[6*pix]; + flush(std::cout); + double xrise[4] = { (cellz[1] - cellz[0])/cell_xlength, // LLL -> HLL + (cellz[3] - cellz[2])/cell_xlength, // LHL -> HHL + (cellz[5] - cellz[4])/cell_xlength, // LLH -> HLH + (cellz[7] - cellz[6])/cell_xlength}; // LHH -> HHH + + double cell_ylength = pillc[6*(pix + numxpill) + 1] - pillc[6*pix + 1]; + double yrise[4] = { (cellz[2] - cellz[0])/cell_ylength, // LLL -> LHL + (cellz[3] - cellz[1])/cell_ylength, // HLL -> HHL + (cellz[6] - cellz[4])/cell_ylength, // LLH -> LHH + (cellz[7] - cellz[5])/cell_ylength}; // HLH -> HHH + + + // Now ignore those edges that touch the global top or bottom surface + // of the entire grdecl model. This is to avoid bias, as these edges probably + // don't follow an overall dip for the model if it exists. + int x_edges = 4; + int y_edges = 4; + std::array gridlimits = getGridLimits(); + double zmin = gridlimits[4]; + double zmax = gridlimits[5]; + // LLL -> HLL + if ((cellz[1] == zmin) || (cellz[0] == zmin)) { + xrise[0] = 0; x_edges--; + } + // LHL -> HHL + if ((cellz[2] == zmin) || (cellz[3] == zmin)) { + xrise[1] = 0; x_edges--; + } + // LLH -> HLH + if ((cellz[4] == zmax) || (cellz[5] == zmax)) { + xrise[2] = 0; x_edges--; + } + // LHH -> HHH + if ((cellz[6] == zmax) || (cellz[7] == zmax)) { + xrise[3] = 0; x_edges--; + } + // LLL -> LHL + if ((cellz[0] == zmin) || (cellz[2] == zmin)) { + yrise[0] = 0; y_edges--; + } + // HLL -> HHL + if ((cellz[1] == zmin) || (cellz[3] == zmin)) { + yrise[1] = 0; y_edges--; + } + // LLH -> LHH + if ((cellz[6] == zmax) || (cellz[4] == zmax)) { + yrise[2] = 0; y_edges--; + } + // HLH -> HHH + if ((cellz[7] == zmax) || (cellz[5] == zmax)) { + yrise[3] = 0; y_edges--; + } + + return std::make_pair( (xrise[0] + xrise[1] + xrise[2] + xrise[3])/x_edges, + (yrise[0] + yrise[1] + yrise[2] + yrise[3])/y_edges); +} +/** + Wrapper for cellDips(i, j, k). +*/ +std::pair EclipseGridInspector::cellDips(int cell_idx) const +{ + std::array idxs = cellIdxToLogicalCoords(cell_idx); + return cellDips(idxs[0], idxs[1], idxs[2]); +} + +std::array EclipseGridInspector::cellIdxToLogicalCoords(int cell_idx) const +{ + + int i,j,k; // Position of cell in cell hierarchy + int horIdx = (cell_idx+1) - int(std::floor(((double)(cell_idx+1))/((double)(logical_gridsize_[0]*logical_gridsize_[1]))))*logical_gridsize_[0]*logical_gridsize_[1]; // index in the corresponding horizon + if (horIdx == 0) { + horIdx = logical_gridsize_[0]*logical_gridsize_[1]; + } + i = horIdx - int(std::floor(((double)horIdx)/((double)logical_gridsize_[0])))*logical_gridsize_[0]; + if (i == 0) { + i = logical_gridsize_[0]; + } + j = (horIdx-i)/logical_gridsize_[0]+1; + k = ((cell_idx+1)-logical_gridsize_[0]*(j-1)-1)/(logical_gridsize_[0]*logical_gridsize_[1])+1; + + std::array a = {{i-1, j-1, k-1}}; + return a; //std::array {{i-1, j-1, k-1}}; +} + +double EclipseGridInspector::cellVolumeVerticalPillars(int i, int j, int k) const +{ + // Checking parameters and obtaining values from parser. + checkLogicalCoords(i, j, k); + const std::vector& pillc = + deck_->getKeyword("COORD").getSIDoubleData(); + int num_pillars = (logical_gridsize_[0] + 1)*(logical_gridsize_[1] + 1); + if (6*num_pillars != int(pillc.size())) { + throw std::runtime_error("Wrong size of COORD field."); + } + const std::vector& z = + deck_->getKeyword("ZCORN").getSIDoubleData(); + int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2]; + if (8*num_cells != int(z.size())) { + throw std::runtime_error("Wrong size of ZCORN field"); + } + + // Computing the base area as half the 2d cross product of the diagonals. + int numxpill = logical_gridsize_[0] + 1; + int pix = i + j*numxpill; + double px[4] = { pillc[6*pix], + pillc[6*(pix + 1)], + pillc[6*(pix + numxpill)], + pillc[6*(pix + numxpill + 1)] }; + double py[4] = { pillc[6*pix + 1], + pillc[6*(pix + 1) + 1], + pillc[6*(pix + numxpill) + 1], + pillc[6*(pix + numxpill + 1) + 1] }; + double diag1[2] = { px[3] - px[0], py[3] - py[0] }; + double diag2[2] = { px[2] - px[1], py[2] - py[1] }; + double area = 0.5*(diag1[0]*diag2[1] - diag1[1]*diag2[0]); + + // Computing the average of the z-differences along each pillar. + int delta[3] = { 1, + 2*logical_gridsize_[0], + 4*logical_gridsize_[0]*logical_gridsize_[1] }; + int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]); + double cellz[8] = { z[ix], z[ix + delta[0]], + z[ix + delta[1]], z[ix + delta[1] + delta[0]], + z[ix + delta[2]], z[ix + delta[2] + delta[0]], + z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] }; + double diffz[4] = { cellz[4] - cellz[0], + cellz[5] - cellz[1], + cellz[6] - cellz[2], + cellz[7] - cellz[3] }; + double averzdiff = 0.25*std::accumulate(diffz, diffz + 4, 0.0); + return averzdiff*area; +} + + +double EclipseGridInspector::cellVolumeVerticalPillars(int cell_idx) const +{ + std::array idxs = cellIdxToLogicalCoords(cell_idx); + return cellVolumeVerticalPillars(idxs[0], idxs[1], idxs[2]); +} + +void EclipseGridInspector::checkLogicalCoords(int i, int j, int k) const +{ + if (i < 0 || i >= logical_gridsize_[0]) + throw std::runtime_error("First coordinate out of bounds"); + if (j < 0 || j >= logical_gridsize_[1]) + throw std::runtime_error("Second coordinate out of bounds"); + if (k < 0 || k >= logical_gridsize_[2]) + throw std::runtime_error("Third coordinate out of bounds"); +} + + +std::array EclipseGridInspector::getGridLimits() const +{ + if (! (deck_->hasKeyword("COORD") && deck_->hasKeyword("ZCORN") && deck_->hasKeyword("SPECGRID")) ) { + throw std::runtime_error("EclipseGridInspector: Grid does not have SPECGRID, COORD, and ZCORN, can't find dimensions."); + } + + std::vector coord = deck_->getKeyword("COORD").getSIDoubleData(); + std::vector zcorn = deck_->getKeyword("ZCORN").getSIDoubleData(); + + double xmin = +DBL_MAX; + double xmax = -DBL_MAX; + double ymin = +DBL_MAX; + double ymax = -DBL_MAX; + + + int pillars = (logical_gridsize_[0]+1) * (logical_gridsize_[1]+1); + + for (int pillarindex = 0; pillarindex < pillars; ++pillarindex) { + if (coord[pillarindex * 6 + 0] > xmax) + xmax = coord[pillarindex * 6 + 0]; + if (coord[pillarindex * 6 + 0] < xmin) + xmin = coord[pillarindex * 6 + 0]; + if (coord[pillarindex * 6 + 1] > ymax) + ymax = coord[pillarindex * 6 + 1]; + if (coord[pillarindex * 6 + 1] < ymin) + ymin = coord[pillarindex * 6 + 1]; + if (coord[pillarindex * 6 + 3] > xmax) + xmax = coord[pillarindex * 6 + 3]; + if (coord[pillarindex * 6 + 3] < xmin) + xmin = coord[pillarindex * 6 + 3]; + if (coord[pillarindex * 6 + 4] > ymax) + ymax = coord[pillarindex * 6 + 4]; + if (coord[pillarindex * 6 + 4] < ymin) + ymin = coord[pillarindex * 6 + 4]; + } + + std::array gridlimits = {{ xmin, xmax, ymin, ymax, + *min_element(zcorn.begin(), zcorn.end()), + *max_element(zcorn.begin(), zcorn.end()) }}; + return gridlimits; +} + + + +std::array EclipseGridInspector::gridSize() const +{ + std::array retval = {{ logical_gridsize_[0], + logical_gridsize_[1], + logical_gridsize_[2] }}; + return retval; +} + + +std::array EclipseGridInspector::cellZvals(int i, int j, int k) const +{ + // Get the zcorn field. + const std::vector& z = deck_->getKeyword("ZCORN").getSIDoubleData(); + int num_cells = logical_gridsize_[0]*logical_gridsize_[1]*logical_gridsize_[2]; + if (8*num_cells != int(z.size())) { + throw std::runtime_error("Wrong size of ZCORN field"); + } + + // Make the coordinate array. + int delta[3] = { 1, + 2*logical_gridsize_[0], + 4*logical_gridsize_[0]*logical_gridsize_[1] }; + int ix = 2*(i*delta[0] + j*delta[1] + k*delta[2]); + std::array cellz = {{ z[ix], z[ix + delta[0]], + z[ix + delta[1]], z[ix + delta[1] + delta[0]], + z[ix + delta[2]], z[ix + delta[2] + delta[0]], + z[ix + delta[2] + delta[1]], z[ix + delta[2] + delta[1] + delta[0]] }}; + return cellz; +} + + +} // namespace Opm diff --git a/src/opm/output/eclipse/EclipseReader.cpp b/src/opm/output/eclipse/EclipseReader.cpp new file mode 100644 index 000000000..6928f9a15 --- /dev/null +++ b/src/opm/output/eclipse/EclipseReader.cpp @@ -0,0 +1,247 @@ +/* + Copyright 2015 Statoil ASA. + + 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 "EclipseReader.hpp" +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include + +#include + +#include + +namespace Opm +{ + + void restoreTemperatureData(const ecl_file_type* file, + EclipseStateConstPtr eclipse_state, + int numcells, + SimulatorState& simulator_state) { + const char* temperature = "TEMP"; + + if (ecl_file_has_kw(file , temperature)) { + ecl_kw_type* temperature_kw = ecl_file_iget_named_kw(file, temperature, 0); + + if (ecl_kw_get_size(temperature_kw) != numcells) { + throw std::runtime_error("Read of restart file: Could not restore temperature data, length of data from file not equal number of cells"); + } + + float* temperature_data = ecl_kw_get_float_ptr(temperature_kw); + + // factor and offset from the temperature values given in the deck to Kelvin + double scaling = eclipse_state->getDeckUnitSystem().parse("Temperature")->getSIScaling(); + double offset = eclipse_state->getDeckUnitSystem().parse("Temperature")->getSIOffset(); + + for (size_t index = 0; index < simulator_state.temperature().size(); ++index) { + simulator_state.temperature()[index] = unit::convert::from((double)temperature_data[index] - offset, scaling); + } + } else { + throw std::runtime_error("Read of restart file: File does not contain TEMP data\n"); + } + } + + + void restorePressureData(const ecl_file_type* file, + EclipseStateConstPtr eclipse_state, + int numcells, + SimulatorState& simulator_state) { + const char* pressure = "PRESSURE"; + + if (ecl_file_has_kw(file , pressure)) { + + ecl_kw_type* pressure_kw = ecl_file_iget_named_kw(file, pressure, 0); + + if (ecl_kw_get_size(pressure_kw) != numcells) { + throw std::runtime_error("Read of restart file: Could not restore pressure data, length of data from file not equal number of cells"); + } + + float* pressure_data = ecl_kw_get_float_ptr(pressure_kw); + const double deck_pressure_unit = (eclipse_state->getDeckUnitSystem().getType() == UnitSystem::UNIT_TYPE_METRIC) ? Opm::unit::barsa : Opm::unit::psia; + for (size_t index = 0; index < simulator_state.pressure().size(); ++index) { + simulator_state.pressure()[index] = unit::convert::from((double)pressure_data[index], deck_pressure_unit); + } + } else { + throw std::runtime_error("Read of restart file: File does not contain PRESSURE data\n"); + } + } + + + void restoreSaturation(const ecl_file_type* file_type, + const PhaseUsage& phaseUsage, + int numcells, + SimulatorState& simulator_state) { + + float* sgas_data = NULL; + float* swat_data = NULL; + + if (phaseUsage.phase_used[BlackoilPhases::Aqua]) { + const char* swat = "SWAT"; + if (ecl_file_has_kw(file_type, swat)) { + ecl_kw_type* swat_kw = ecl_file_iget_named_kw(file_type , swat, 0); + swat_data = ecl_kw_get_float_ptr(swat_kw); + std::vector swat_datavec(&swat_data[0], &swat_data[numcells]); + EclipseIOUtil::addToStripedData(swat_datavec, simulator_state.saturation(), phaseUsage.phase_pos[BlackoilPhases::Aqua], phaseUsage.num_phases); + } else { + std::string error_str = "Restart file is missing SWAT data!\n"; + throw std::runtime_error(error_str); + } + } + + if (phaseUsage.phase_used[BlackoilPhases::Vapour]) { + const char* sgas = "SGAS"; + if (ecl_file_has_kw(file_type, sgas)) { + ecl_kw_type* sgas_kw = ecl_file_iget_named_kw(file_type , sgas, 0); + sgas_data = ecl_kw_get_float_ptr(sgas_kw); + std::vector sgas_datavec(&sgas_data[0], &sgas_data[numcells]); + EclipseIOUtil::addToStripedData(sgas_datavec, simulator_state.saturation(), phaseUsage.phase_pos[BlackoilPhases::Vapour], phaseUsage.num_phases); + } else { + std::string error_str = "Restart file is missing SGAS data!\n"; + throw std::runtime_error(error_str); + } + } + } + + + void restoreRSandRV(const ecl_file_type* file_type, + SimulationConfigConstPtr sim_config, + int numcells, + BlackoilState* blackoil_state) { + + if (sim_config->hasDISGAS()) { + const char* RS = "RS"; + if (ecl_file_has_kw(file_type, RS)) { + ecl_kw_type* rs_kw = ecl_file_iget_named_kw(file_type, RS, 0); + float* rs_data = ecl_kw_get_float_ptr(rs_kw); + std::vector rs_datavec(&rs_data[0], &rs_data[numcells]); + blackoil_state->gasoilratio().clear(); + blackoil_state->gasoilratio().insert(blackoil_state->gasoilratio().begin(), rs_datavec.begin(), rs_datavec.end()); + } else { + throw std::runtime_error("Restart file is missing RS data!\n"); + } + } + + if (sim_config->hasVAPOIL()) { + const char* RV = "RV"; + if (ecl_file_has_kw(file_type, RV)) { + ecl_kw_type* rv_kw = ecl_file_iget_named_kw(file_type, RV, 0); + float* rv_data = ecl_kw_get_float_ptr(rv_kw); + std::vector rv_datavec(&rv_data[0], &rv_data[numcells]); + blackoil_state->rv().clear(); + blackoil_state->rv().insert(blackoil_state->rv().begin(), rv_datavec.begin(), rv_datavec.end()); + } else { + throw std::runtime_error("Restart file is missing RV data!\n"); + } + } + } + + + void restoreSOLUTION(const std::string& restart_filename, + int reportstep, + bool unified, + EclipseStateConstPtr eclipseState, + int numcells, + const PhaseUsage& phaseUsage, + SimulatorState& simulator_state) + { + const char* filename = restart_filename.c_str(); + ecl_file_type* file_type = ecl_file_open(filename, 0); + + if (file_type) { + bool block_selected = unified ? ecl_file_select_rstblock_report_step(file_type , reportstep) : true; + + if (block_selected) { + restorePressureData(file_type, eclipseState, numcells, simulator_state); + restoreTemperatureData(file_type, eclipseState, numcells, simulator_state); + restoreSaturation(file_type, phaseUsage, numcells, simulator_state); + BlackoilState* blackoilState = dynamic_cast(&simulator_state); + if (blackoilState) { + SimulationConfigConstPtr sim_config = eclipseState->getSimulationConfig(); + restoreRSandRV(file_type, sim_config, numcells, blackoilState); + } + } else { + std::string error_str = "Restart file " + restart_filename + " does not contain data for report step " + std::to_string(reportstep) + "!\n"; + throw std::runtime_error(error_str); + } + ecl_file_close(file_type); + } else { + std::string error_str = "Restart file " + restart_filename + " not found!\n"; + throw std::runtime_error(error_str); + } + } + + + void restoreOPM_XWELKeyword(const std::string& restart_filename, int reportstep, bool unified, WellState& wellstate) + { + const char * keyword = "OPM_XWEL"; + const char* filename = restart_filename.c_str(); + ecl_file_type* file_type = ecl_file_open(filename, 0); + + if (file_type != NULL) { + + bool block_selected = unified ? ecl_file_select_rstblock_report_step(file_type , reportstep) : true; + + if (block_selected) { + ecl_kw_type* xwel = ecl_file_iget_named_kw(file_type , keyword, 0); + const double* xwel_data = ecl_kw_get_double_ptr(xwel); + std::copy_n(xwel_data + wellstate.getRestartTemperatureOffset(), wellstate.temperature().size(), wellstate.temperature().begin()); + std::copy_n(xwel_data + wellstate.getRestartBhpOffset(), wellstate.bhp().size(), wellstate.bhp().begin()); + std::copy_n(xwel_data + wellstate.getRestartPerfPressOffset(), wellstate.perfPress().size(), wellstate.perfPress().begin()); + std::copy_n(xwel_data + wellstate.getRestartPerfRatesOffset(), wellstate.perfRates().size(), wellstate.perfRates().begin()); + std::copy_n(xwel_data + wellstate.getRestartWellRatesOffset(), wellstate.wellRates().size(), wellstate.wellRates().begin()); + } else { + std::string error_str = "Restart file " + restart_filename + " does not contain data for report step " + std::to_string(reportstep) + "!\n"; + throw std::runtime_error(error_str); + } + ecl_file_close(file_type); + } else { + std::string error_str = "Restart file " + restart_filename + " not found!\n"; + throw std::runtime_error(error_str); + } + } + + + + void init_from_restart_file(EclipseStateConstPtr eclipse_state, + int numcells, + const PhaseUsage& phase_usage, + SimulatorState& simulator_state, + WellState& wellstate) { + + InitConfigConstPtr initConfig = eclipse_state->getInitConfig(); + IOConfigConstPtr ioConfig = eclipse_state->getIOConfig(); + int restart_step = initConfig->getRestartStep(); + const std::string& restart_file_root = initConfig->getRestartRootName(); + bool output = false; + const std::string& restart_file_name = ioConfig->getRestartFileName(restart_file_root, restart_step, output); + + Opm::restoreSOLUTION(restart_file_name, restart_step, ioConfig->getUNIFIN(), eclipse_state, numcells, phase_usage, simulator_state); + Opm::restoreOPM_XWELKeyword(restart_file_name, restart_step, ioConfig->getUNIFIN(), wellstate); + } + + +} // namespace Opm diff --git a/src/opm/output/eclipse/EclipseWriteRFTHandler.cpp b/src/opm/output/eclipse/EclipseWriteRFTHandler.cpp new file mode 100644 index 000000000..17fb89030 --- /dev/null +++ b/src/opm/output/eclipse/EclipseWriteRFTHandler.cpp @@ -0,0 +1,146 @@ +/* + Copyright 2015 Statoil ASA. + + 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 + + + +namespace Opm { +namespace EclipseWriterDetails { + + EclipseWriteRFTHandler::EclipseWriteRFTHandler(const int * compressedToCartesianCellIdx, size_t numCells, size_t cartesianSize) { + initGlobalToActiveIndex(compressedToCartesianCellIdx, numCells, cartesianSize); + } + + void EclipseWriteRFTHandler::writeTimeStep(const std::string& filename, + const ert_ecl_unit_enum ecl_unit, + const SimulatorTimerInterface& simulatorTimer, + std::vector& wells, + EclipseGridConstPtr eclipseGrid, + std::vector& pressure, + std::vector& swat, + std::vector& sgas) { + + + + std::vector rft_nodes; + for (std::vector::const_iterator ci = wells.begin(); ci != wells.end(); ++ci) { + WellConstPtr well = *ci; + if ((well->getRFTActive(simulatorTimer.reportStepNum())) || (well->getPLTActive(simulatorTimer.reportStepNum()))) { + ecl_rft_node_type * ecl_node = createEclRFTNode(well, + simulatorTimer, + eclipseGrid, + pressure, + swat, + sgas); + + // TODO: replace this silenced warning with an appropriate + // use of the OpmLog facilities. + // if (well->getPLTActive(simulatorTimer.reportStepNum())) { + // std::cerr << "PLT not supported, writing RFT data" << std::endl; + // } + + rft_nodes.push_back(ecl_node); + } + } + + + if (rft_nodes.size() > 0) { + ecl_rft_file_update(filename.c_str(), rft_nodes.data(), rft_nodes.size(), ecl_unit); + } + + //Cleanup: The ecl_rft_file_update method takes care of freeing the ecl_rft_nodes that it receives. + // Each ecl_rft_node is again responsible for freeing it's cells. + } + + + + + ecl_rft_node_type * EclipseWriteRFTHandler::createEclRFTNode(WellConstPtr well, + const SimulatorTimerInterface& simulatorTimer, + EclipseGridConstPtr eclipseGrid, + const std::vector& pressure, + const std::vector& swat, + const std::vector& sgas) { + + + const std::string& well_name = well->name(); + size_t timestep = (size_t)simulatorTimer.reportStepNum(); + time_t recording_date = simulatorTimer.currentPosixTime(); + double days = Opm::unit::convert::to(simulatorTimer.simulationTimeElapsed(), Opm::unit::day); + + std::string type = "RFT"; + ecl_rft_node_type * ecl_rft_node = ecl_rft_node_alloc_new(well_name.c_str(), type.c_str(), recording_date, days); + + CompletionSetConstPtr completionsSet = well->getCompletions(timestep); + for (size_t index = 0; index < completionsSet->size(); ++index) { + CompletionConstPtr completion = completionsSet->get(index); + size_t i = (size_t)completion->getI(); + size_t j = (size_t)completion->getJ(); + size_t k = (size_t)completion->getK(); + + size_t global_index = eclipseGrid->getGlobalIndex(i,j,k); + int active_index = globalToActiveIndex_[global_index]; + + if (active_index > -1) { + double depth = eclipseGrid->getCellDepth(i,j,k); + double completion_pressure = pressure.size() > 0 ? pressure[active_index] : 0.0; + double saturation_water = swat.size() > 0 ? swat[active_index] : 0.0; + double saturation_gas = sgas.size() > 0 ? sgas[active_index] : 0.0; + + ecl_rft_cell_type * ecl_rft_cell = ecl_rft_cell_alloc_RFT( i ,j, k , depth, completion_pressure, saturation_water, saturation_gas); + ecl_rft_node_append_cell( ecl_rft_node , ecl_rft_cell); + } + } + + return ecl_rft_node; + } + + + void EclipseWriteRFTHandler::initGlobalToActiveIndex(const int * compressedToCartesianCellIdx, size_t numCells, size_t cartesianSize) { + globalToActiveIndex_.resize(cartesianSize, -1); + for (size_t active_index = 0; active_index < numCells; ++active_index) { + //If compressedToCartesianCellIdx is NULL, assume no compressed to cartesian mapping, set global equal to active index + int global_index = (NULL != compressedToCartesianCellIdx) ? compressedToCartesianCellIdx[active_index] : active_index; + globalToActiveIndex_[global_index] = active_index; + } + } + + +}//namespace EclipseWriterDetails +}//namespace Opm diff --git a/src/opm/output/eclipse/EclipseWriter.cpp b/src/opm/output/eclipse/EclipseWriter.cpp new file mode 100644 index 000000000..98ecf42ea --- /dev/null +++ b/src/opm/output/eclipse/EclipseWriter.cpp @@ -0,0 +1,1501 @@ +/* + 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 "config.h" + +#include "EclipseWriter.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include // WellType + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include // to_upper_copy +#include +#include // path + +#include // mktime +#include +#include // unique_ptr +#include // move + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#define OPM_XWEL "OPM_XWEL" + +// namespace start here since we don't want the ERT headers in it +namespace Opm { +namespace EclipseWriterDetails { +/// Names of the saturation property for each phase. The order of these +/// names are critical; they must be the same as the BlackoilPhases enum +static const char* saturationKeywordNames[] = { "SWAT", "SOIL", "SGAS" }; + +// throw away the data for all non-active cells and reorder to the Cartesian logic of +// eclipse +void restrictAndReorderToActiveCells(std::vector &data, + int numCells, + const int* compressedToCartesianCellIdx) +{ + if (!compressedToCartesianCellIdx) + // if there is no active -> global mapping, all cells + // are considered active + return; + + std::vector eclData; + eclData.reserve( numCells ); + + // activate those cells that are actually there + for (int i = 0; i < numCells; ++i) { + eclData.push_back( data[ compressedToCartesianCellIdx[i] ] ); + } + data.swap( eclData ); +} + +// convert the units of an array +void convertFromSiTo(std::vector &siValues, double toSiConversionFactor, double toSiOffset = 0) +{ + for (size_t curIdx = 0; curIdx < siValues.size(); ++curIdx) { + siValues[curIdx] = unit::convert::to(siValues[curIdx] - toSiOffset, toSiConversionFactor); + } +} + +// extract a sub-array of a larger one which represents multiple +// striped ones +void extractFromStripedData(std::vector &data, + int offset, + int stride) +{ + size_t tmpIdx = 0; + for (size_t curIdx = offset; curIdx < data.size(); curIdx += stride) { + assert(tmpIdx <= curIdx); + data[tmpIdx] = data[curIdx]; + ++tmpIdx; + } + // shirk the result + data.resize(tmpIdx); +} + +/// Convert OPM phase usage to ERT bitmask +int ertPhaseMask(const PhaseUsage uses) +{ + return (uses.phase_used[BlackoilPhases::Liquid] ? ECL_OIL_PHASE : 0) + | (uses.phase_used[BlackoilPhases::Aqua] ? ECL_WATER_PHASE : 0) + | (uses.phase_used[BlackoilPhases::Vapour] ? ECL_GAS_PHASE : 0); +} + + +/** + * Eclipse "keyword" (i.e. named data) for a vector. + */ +template +class Keyword : private boost::noncopyable +{ +public: + // Default constructor + Keyword() + : ertHandle_(0) + {} + + /// Initialization from double-precision array. + Keyword(const std::string& name, + const std::vector& data) + : ertHandle_(0) + { set(name, data); } + + /// Initialization from double-precision array. + Keyword(const std::string& name, + const std::vector& data) + : ertHandle_(0) + { set(name, data); } + + Keyword(const std::string& name, + const std::vector& data) + : ertHandle_(0) + {set(name, data); } + + ~Keyword() + { + if (ertHandle_) + ecl_kw_free(ertHandle_); + } + + template + void set(const std::string name, const std::vector& data) + { + + if(ertHandle_) { + ecl_kw_free(ertHandle_); + } + + + ertHandle_ = ecl_kw_alloc(name.c_str(), + data.size(), + ertType_()); + + // number of elements to take + const int numEntries = data.size(); + + // fill it with values + + T* target = static_cast(ecl_kw_get_ptr(ertHandle())); + for (int i = 0; i < numEntries; ++i) { + target[i] = static_cast(data[i]); + } + } + + void set(const std::string name, const std::vector& data) + { + if(ertHandle_) { + ecl_kw_free(ertHandle_); + } + + + ertHandle_ = ecl_kw_alloc(name.c_str(), + data.size(), + ertType_()); + + // number of elements to take + const int numEntries = data.size(); + for (int i = 0; i < numEntries; ++i) { + ecl_kw_iset_char_ptr( ertHandle_, i, data[i]); + } + + } + + ecl_kw_type *ertHandle() const + { return ertHandle_; } + +private: + static ecl_type_enum ertType_() + { + if (std::is_same::value) + { return ECL_FLOAT_TYPE; } + if (std::is_same::value) + { return ECL_DOUBLE_TYPE; } + if (std::is_same::value) + { return ECL_INT_TYPE; } + if (std::is_same::value) + { return ECL_CHAR_TYPE; } + + + OPM_THROW(std::logic_error, + "Unhandled type for data elements in EclipseWriterDetails::Keyword"); + } + + ecl_kw_type *ertHandle_; +}; + +/** + * Pointer to memory that holds the name to an Eclipse output file. + */ +class FileName : private boost::noncopyable +{ +public: + FileName(const std::string& outputDir, + const std::string& baseName, + ecl_file_enum type, + int writeStepIdx, + bool formatted) + { + ertHandle_ = ecl_util_alloc_filename(outputDir.c_str(), + baseName.c_str(), + type, + formatted, + writeStepIdx); + } + + ~FileName() + { std::free(ertHandle_); } + + const char *ertHandle() const + { return ertHandle_; } + +private: + char *ertHandle_; +}; + +class Restart : private boost::noncopyable +{ +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, + IOConfigConstPtr ioConfig) + { + + ecl_file_enum type_of_restart_file = ioConfig->getUNIFOUT() ? ECL_UNIFIED_RESTART_FILE : ECL_RESTART_FILE; + + restartFileName_ = ecl_util_alloc_filename(outputDir.c_str(), + baseName.c_str(), + type_of_restart_file, + ioConfig->getFMTOUT(), // use formatted instead of binary output? + writeStepIdx); + + if ((writeStepIdx > 0) && (ECL_UNIFIED_RESTART_FILE == type_of_restart_file)) { + restartFileHandle_ = ecl_rst_file_open_append(restartFileName_); + } + else { + restartFileHandle_ = ecl_rst_file_open_write(restartFileName_); + } + } + + template + void add_kw(const Keyword& kw) + { ecl_rst_file_add_kw(restartFileHandle_, kw.ertHandle()); } + + + void addRestartFileIwelData(std::vector& iwel_data, size_t currentStep , WellConstPtr well, size_t offset) const { + CompletionSetConstPtr completions = well->getCompletions( currentStep ); + + iwel_data[ offset + IWEL_HEADI_ITEM ] = well->getHeadI() + 1; + iwel_data[ offset + IWEL_HEADJ_ITEM ] = well->getHeadJ() + 1; + iwel_data[ offset + IWEL_CONNECTIONS_ITEM ] = completions->size(); + iwel_data[ offset + IWEL_GROUP_ITEM ] = 1; + + { + WellType welltype = well->isProducer(currentStep) ? PRODUCER : INJECTOR; + int ert_welltype = EclipseWriter::eclipseWellTypeMask(welltype, well->getInjectionProperties(currentStep).injectorType); + iwel_data[ offset + IWEL_TYPE_ITEM ] = ert_welltype; + } + + iwel_data[ offset + IWEL_STATUS_ITEM ] = EclipseWriter::eclipseWellStatusMask(well->getStatus(currentStep)); + } + + void addRestartFileXwelData(const WellState& wellState, std::vector& xwel_data) const { + std::copy_n(wellState.bhp().begin(), wellState.bhp().size(), xwel_data.begin() + wellState.getRestartBhpOffset()); + std::copy_n(wellState.perfPress().begin(), wellState.perfPress().size(), xwel_data.begin() + wellState.getRestartPerfPressOffset()); + std::copy_n(wellState.perfRates().begin(), wellState.perfRates().size(), xwel_data.begin() + wellState.getRestartPerfRatesOffset()); + std::copy_n(wellState.temperature().begin(), wellState.temperature().size(), xwel_data.begin() + wellState.getRestartTemperatureOffset()); + std::copy_n(wellState.wellRates().begin(), wellState.wellRates().size(), xwel_data.begin() + wellState.getRestartWellRatesOffset()); + } + + void addRestartFileIconData(std::vector& icon_data, CompletionSetConstPtr completions , size_t wellICONOffset) const { + for (size_t i = 0; i < completions->size(); ++i) { + CompletionConstPtr completion = completions->get(i); + size_t iconOffset = wellICONOffset + i * Opm::EclipseWriterDetails::Restart::NICONZ; + icon_data[ iconOffset + ICON_IC_ITEM] = 1; + + icon_data[ iconOffset + ICON_I_ITEM] = completion->getI() + 1; + icon_data[ iconOffset + ICON_J_ITEM] = completion->getJ() + 1; + icon_data[ iconOffset + ICON_K_ITEM] = completion->getK() + 1; + + { + WellCompletion::StateEnum completion_state = completion->getState(); + if (completion_state == WellCompletion::StateEnum::OPEN) { + icon_data[ iconOffset + ICON_STATUS_ITEM ] = 1; + } else { + icon_data[ iconOffset + ICON_STATUS_ITEM ] = 0; + } + } + icon_data[ iconOffset + ICON_DIRECTION_ITEM] = (int)completion->getDirection(); + } + } + + + ~Restart() + { + free(restartFileName_); + ecl_rst_file_close(restartFileHandle_); + } + + void writeHeader(const SimulatorTimerInterface& /*timer*/, + int writeStepIdx, + 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(restartFileHandle_, + writeStepIdx, + rsthead_data); + + } + + ecl_rst_file_type *ertHandle() const + { return restartFileHandle_; } + +private: + char *restartFileName_; + ecl_rst_file_type *restartFileHandle_; +}; + +/** + * 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 : private boost::noncopyable +{ +public: + Solution(Restart& restartHandle) + : restartHandle_(&restartHandle) + { ecl_rst_file_start_solution(restartHandle_->ertHandle()); } + + ~Solution() + { ecl_rst_file_end_solution(restartHandle_->ertHandle()); } + + template + void add(const Keyword& kw) + { ecl_rst_file_add_kw(restartHandle_->ertHandle(), kw.ertHandle()); } + + ecl_rst_file_type *ertHandle() const + { return restartHandle_->ertHandle(); } + +private: + Restart* restartHandle_; +}; + +/// Supported well types. Enumeration doesn't let us get all the members, +/// so we must have an explicit array. +static WellType WELL_TYPES[] = { INJECTOR, PRODUCER }; + +class WellReport; + +class Summary : private boost::noncopyable +{ +public: + Summary(const std::string& outputDir, + const std::string& baseName, + const SimulatorTimerInterface& timer, + bool time_in_days , + int nx, + int ny, + int nz) + { + boost::filesystem::path casePath(outputDir); + casePath /= boost::to_upper_copy(baseName); + + // convert the start time to seconds since 1970-1-1@00:00:00 + boost::posix_time::ptime startTime + = timer.startDateTime(); + tm t = boost::posix_time::to_tm(startTime); + double secondsSinceEpochStart = std::mktime(&t); + + ertHandle_ = ecl_sum_alloc_writer(casePath.string().c_str(), + false, /* formatted */ + true, /* unified */ + ":", /* join string */ + secondsSinceEpochStart, + time_in_days, + nx, + ny, + nz); + } + + ~Summary() + { ecl_sum_free(ertHandle_); } + + typedef std::unique_ptr SummaryReportVar; + typedef std::vector SummaryReportVarCollection; + + Summary& addWell(SummaryReportVar var) + { + summaryReportVars_.push_back(std::move(var)); + return *this; + } + + // no inline implementation of these two methods since they depend + // on the classes defined in the following. + + // add rate variables for each of the well in the input file + void addAllWells(Opm::EclipseStateConstPtr eclipseState, + const PhaseUsage& uses); + void writeTimeStep(int writeStepIdx, + const SimulatorTimerInterface& timer, + const WellState& wellState); + + ecl_sum_type *ertHandle() const + { return ertHandle_; } + +private: + ecl_sum_type *ertHandle_; + + Opm::EclipseStateConstPtr eclipseState_; + SummaryReportVarCollection summaryReportVars_; +}; + +class SummaryTimeStep : private boost::noncopyable +{ +public: + SummaryTimeStep(Summary& summaryHandle, + int writeStepIdx, + const SimulatorTimerInterface &timer) + { + ertHandle_ = ecl_sum_add_tstep(summaryHandle.ertHandle(), + writeStepIdx, + timer.simulationTimeElapsed()); + } + + // no destructor in this class as ERT takes care of freeing the + // handle as part of freeing the solution handle! + + ecl_sum_tstep_type *ertHandle() const + { return ertHandle_; }; + +private: + ecl_sum_tstep_type *ertHandle_; +}; + + +/** + * Initialization file which contains static properties (such as + * porosity and permeability) for the simulation field. + */ +class Init : private boost::noncopyable +{ +public: + Init(const std::string& outputDir, + const std::string& baseName, + int writeStepIdx, + IOConfigConstPtr ioConfig) : egridFileName_(outputDir, + baseName, + ECL_EGRID_FILE, + writeStepIdx, + ioConfig->getFMTOUT()) + { + bool formatted = ioConfig->getFMTOUT(); + + FileName initFileName(outputDir, + baseName, + ECL_INIT_FILE, + writeStepIdx, + formatted); + + ertHandle_ = fortio_open_writer(initFileName.ertHandle(), + formatted, + ECL_ENDIAN_FLIP); + } + + ~Init() + { fortio_fclose(ertHandle_); } + + void writeHeader(int numCells, + const int* compressedToCartesianCellIdx, + const SimulatorTimerInterface& timer, + Opm::EclipseStateConstPtr eclipseState, + const PhaseUsage uses) + { + auto dataField = eclipseState->getDoubleGridProperty("PORO")->getData(); + restrictAndReorderToActiveCells(dataField, numCells, compressedToCartesianCellIdx); + + auto eclGrid = eclipseState->getEclipseGridCopy(); + + // update the ACTNUM array using the processed cornerpoint grid + std::vector actnumData(eclGrid->getCartesianSize(), 1); + if (compressedToCartesianCellIdx) { + std::fill(actnumData.begin(), actnumData.end(), 0); + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + int cartesianCellIdx = compressedToCartesianCellIdx[cellIdx]; + actnumData[cartesianCellIdx] = 1; + } + } + + eclGrid->resetACTNUM(&actnumData[0]); + + + // finally, write the grid to disk + IOConfigConstPtr ioConfig = eclipseState->getIOConfigConst(); + if (ioConfig->getWriteEGRIDFile()) { + if (eclipseState->getDeckUnitSystem().getType() == UnitSystem::UNIT_TYPE_METRIC){ + eclGrid->fwriteEGRID(egridFileName_.ertHandle(), true); + }else{ + eclGrid->fwriteEGRID(egridFileName_.ertHandle(), false); + } + } + + + if (ioConfig->getWriteINITFile()) { + Keyword poro_kw("PORO", dataField); + ecl_init_file_fwrite_header(ertHandle(), + eclGrid->c_ptr(), + poro_kw.ertHandle(), + ertPhaseMask(uses), + timer.currentPosixTime()); + } + } + + void writeKeyword(const std::string& keywordName, const std::vector &data) + { + Keyword kw(keywordName, data); + ecl_kw_fwrite(kw.ertHandle(), ertHandle()); + } + + fortio_type *ertHandle() const + { return ertHandle_; } + +private: + fortio_type *ertHandle_; + FileName egridFileName_; +}; + + + + +/** + * Summary variable that reports a characteristics of a well. + */ +class WellReport : private boost::noncopyable +{ +protected: + WellReport(const Summary& summary, /* section to add to */ + Opm::EclipseStateConstPtr eclipseState, + Opm::WellConstPtr& well, + PhaseUsage uses, /* phases present */ + BlackoilPhases::PhaseIndex phase, /* oil, water or gas */ + WellType type, /* prod. or inj. */ + char aggregation, /* rate or total */ + std::string unit) + // save these for when we update the value in a timestep + : eclipseState_(eclipseState) + , well_(well) + , phaseUses_(uses) + , phaseIdx_(phase) + { + // producers can be seen as negative injectors + if (type == INJECTOR) + sign_ = +1.0; + else + sign_ = -1.0; + ertHandle_ = ecl_sum_add_var(summary.ertHandle(), + varName_(phase, + type, + aggregation).c_str(), + well_->name().c_str(), + /*num=*/ 0, + unit.c_str(), + /*defaultValue=*/ 0.); + } + +public: + /// Retrieve the value which the monitor is supposed to write to the summary file + /// according to the state of the well. + virtual double retrieveValue(const int writeStepIdx, + const SimulatorTimerInterface& timer, + const WellState& wellState, + const std::map& nameToIdxMap) = 0; + + smspec_node_type *ertHandle() const + { return ertHandle_; } + +protected: + + + void updateTimeStepWellIndex_(const std::map& nameToIdxMap) + { + const std::string& wellName = well_->name(); + + const auto wellIdxIt = nameToIdxMap.find(wellName); + if (wellIdxIt == nameToIdxMap.end()) { + timeStepWellIdx_ = -1; + flatIdx_ = -1; + return; + } + + timeStepWellIdx_ = wellIdxIt->second; + flatIdx_ = timeStepWellIdx_*phaseUses_.num_phases + phaseUses_.phase_pos[phaseIdx_]; + } + + // return m^3/s of injected or produced fluid + double rate(const WellState& wellState) + { + double value = 0; + if (wellState.wellRates().size() > 0) { + assert(int(wellState.wellRates().size()) > flatIdx_); + value = sign_ * wellState.wellRates()[flatIdx_]; + } + return value; + } + + double bhp(const WellState& wellState) + { + if (wellState.bhp().size() > 0) { + // Note that 'flatIdx_' is used here even though it is meant + // to give a (well,phase) pair. + const int numPhases = wellState.wellRates().size() / wellState.bhp().size(); + + return wellState.bhp()[flatIdx_/numPhases]; + } + return 0.0; + } + + /// Get the index associated a well name + int wellIndex_(Opm::EclipseStateConstPtr eclipseState) + { + const Opm::ScheduleConstPtr schedule = eclipseState->getSchedule(); + + const std::string& wellName = well_->name(); + const auto& wells = schedule->getWells(); + for (size_t wellIdx = 0; wellIdx < wells.size(); ++wellIdx) { + if (wells[wellIdx]->name() == wellName) { + return wellIdx; + } + } + + OPM_THROW(std::runtime_error, + "Well '" << wellName << "' is not present in deck"); + } + + /// Compose the name of the summary variable, e.g. "WOPR" for + /// well oil production rate. + std::string varName_(BlackoilPhases::PhaseIndex phase, + WellType type, + char aggregation) + { + std::string name; + name += 'W'; // well + if (aggregation == 'B') { + name += "BHP"; + } else { + switch (phase) { + case BlackoilPhases::Aqua: name += 'W'; break; /* water */ + case BlackoilPhases::Vapour: name += 'G'; break; /* gas */ + case BlackoilPhases::Liquid: name += 'O'; break; /* oil */ + default: + OPM_THROW(std::runtime_error, + "Unknown phase used in blackoil reporting"); + } + switch (type) { + case WellType::INJECTOR: name += 'I'; break; + case WellType::PRODUCER: name += 'P'; break; + default: + OPM_THROW(std::runtime_error, + "Unknown well type used in blackoil reporting"); + } + name += aggregation; /* rate ('R') or total ('T') */ + } + return name; + } + + smspec_node_type *ertHandle_; + + Opm::EclipseStateConstPtr eclipseState_; + Opm::WellConstPtr well_; + + PhaseUsage phaseUses_; + BlackoilPhases::PhaseIndex phaseIdx_; + + int timeStepWellIdx_; + + /// index into a (flattened) wellsOfTimeStep*phases matrix + int flatIdx_; + + /// natural sign of the rate + double sign_; +}; + +/// Monitors the rate given by a well. +class WellRate : public WellReport +{ +public: + WellRate(const Summary& summary, + Opm::EclipseStateConstPtr eclipseState, + Opm::WellConstPtr well, + PhaseUsage uses, + BlackoilPhases::PhaseIndex phase, + WellType type, + UnitSystem::UnitType unitType) + : WellReport(summary, + eclipseState, + well, + uses, + phase, + type, + 'R', + handleUnit_(phase, unitType)) + { + } + + virtual double retrieveValue(const int /* writeStepIdx */, + const SimulatorTimerInterface& timer, + const WellState& wellState, + const std::map& wellNameToIdxMap) + { + // find the index for the quantity in the wellState + this->updateTimeStepWellIndex_(wellNameToIdxMap); + if (this->flatIdx_ < 0) { + // well not active in current time step + return 0.0; + } + + if (well_->getStatus(timer.reportStepNum()) == WellCommon::SHUT) { + // well is shut in the current time step + return 0.0; + } + + // TODO: Why only positive rates? + using namespace Opm::unit; + return convert::to(std::max(0., rate(wellState)), + targetRateToSiConversionFactor_); + } + +private: + const std::string handleUnit_(BlackoilPhases::PhaseIndex phase, UnitSystem::UnitType unitType) { + using namespace Opm::unit; + if (phase == BlackoilPhases::Liquid || phase == BlackoilPhases::Aqua) { + if (unitType == UnitSystem::UNIT_TYPE_FIELD) { + unitName_ = "STB/DAY"; + targetRateToSiConversionFactor_ = stb/day; // m^3/s -> STB/day + } + else if (unitType == UnitSystem::UNIT_TYPE_METRIC) { + unitName_ = "SM3/DAY"; + targetRateToSiConversionFactor_ = cubic(meter)/day; // m^3/s -> m^3/day + } + else + OPM_THROW(std::logic_error, "Deck uses unexpected unit system"); + } + else if (phase == BlackoilPhases::Vapour) { + if (unitType == UnitSystem::UNIT_TYPE_FIELD) { + unitName_ = "MSCF/DAY"; + targetRateToSiConversionFactor_ = 1000*cubic(feet)/day; // m^3/s -> MSCF^3/day + } + else if (unitType == UnitSystem::UNIT_TYPE_METRIC) { + unitName_ = "SM3/DAY"; + targetRateToSiConversionFactor_ = cubic(meter)/day; // m^3/s -> m^3/day + } + else + OPM_THROW(std::logic_error, "Deck uses unexpected unit system"); + } + else + OPM_THROW(std::logic_error, + "Unexpected phase " << phase); + return unitName_; + } + + const char* unitName_; + double targetRateToSiConversionFactor_; +}; + +/// Monitors the total production in a well. +class WellTotal : public WellReport +{ +public: + WellTotal(const Summary& summary, + Opm::EclipseStateConstPtr eclipseState, + Opm::WellConstPtr well, + PhaseUsage uses, + BlackoilPhases::PhaseIndex phase, + WellType type, + UnitSystem::UnitType unitType) + : WellReport(summary, + eclipseState, + well, + uses, + phase, + type, + 'T', + handleUnit_(phase, unitType)) + // nothing produced when the reporting starts + , total_(0.) + { } + + virtual double retrieveValue(const int writeStepIdx, + const SimulatorTimerInterface& timer, + const WellState& wellState, + const std::map& wellNameToIdxMap) + { + if (writeStepIdx == 0) { + // We are at the initial state. + // No step has been taken yet. + return 0.0; + } + + if (well_->getStatus(timer.reportStepNum()) == WellCommon::SHUT) { + // well is shut in the current time step + return 0.0; + } + + // find the index for the quantity in the wellState + this->updateTimeStepWellIndex_(wellNameToIdxMap); + if (this->flatIdx_ < 0) { + // well not active in current time step + return 0.0; + } + + // due to using an Euler method as time integration scheme, the well rate is the + // average for the time step. For more complicated time stepping schemes, the + // integral of the rate is not simply multiplying two numbers... + const double intg = timer.stepLengthTaken() * rate(wellState); + + // add this timesteps production to the total + total_ += intg; + // report the new production total + return unit::convert::to(total_, targetRateToSiConversionFactor_); + } + +private: + const std::string handleUnit_(BlackoilPhases::PhaseIndex phase, UnitSystem::UnitType unitType) { + using namespace Opm::unit; + if (phase == BlackoilPhases::Liquid || phase == BlackoilPhases::Aqua) { + if (unitType == UnitSystem::UNIT_TYPE_FIELD) { + unitName_ = "STB"; + targetRateToSiConversionFactor_ = stb; // m^3 -> STB + } + else if (unitType == UnitSystem::UNIT_TYPE_METRIC) { + unitName_ = "SM3"; + targetRateToSiConversionFactor_ = cubic(meter); // m^3 -> m^3 + } + else + OPM_THROW(std::logic_error, "Deck uses unexpected unit system"); + } + else if (phase == BlackoilPhases::Vapour) { + if (unitType == UnitSystem::UNIT_TYPE_FIELD) { + unitName_ = "MSCF"; + targetRateToSiConversionFactor_ = 1000*cubic(feet); // m^3 -> MSCF^3 + } + else if (unitType == UnitSystem::UNIT_TYPE_METRIC) { + unitName_ = "SM3"; + targetRateToSiConversionFactor_ = cubic(meter); // m^3 -> m^3 + } + else + OPM_THROW(std::logic_error, "Deck uses unexpected unit system"); + } + else + OPM_THROW(std::logic_error, + "Unexpected phase " << phase); + return unitName_; + } + + const char* unitName_; + double targetRateToSiConversionFactor_; + + /// Aggregated value of the course of the simulation + double total_; +}; + +/// Monitors the bottom hole pressure in a well. +class WellBhp : public WellReport +{ +public: + WellBhp(const Summary& summary, + Opm::EclipseStateConstPtr eclipseState, + Opm::WellConstPtr well, + PhaseUsage uses, + BlackoilPhases::PhaseIndex phase, + WellType type, + UnitSystem::UnitType unitType) + : WellReport(summary, + eclipseState, + well, + uses, + phase, + type, + 'B', + handleUnit_(unitType)) + { } + + virtual double retrieveValue(const int /* writeStepIdx */, + const SimulatorTimerInterface& timer, + const WellState& wellState, + const std::map& wellNameToIdxMap) + { + // find the index for the quantity in the wellState + this->updateTimeStepWellIndex_(wellNameToIdxMap); + if (this->flatIdx_ < 0) { + // well not active in current time step + return 0.0; + } + if (well_->getStatus(timer.reportStepNum()) == WellCommon::SHUT) { + // well is shut in the current time step + return 0.0; + } + + return unit::convert::to(bhp(wellState), targetRateToSiConversionFactor_); + } + +private: + const std::string handleUnit_(UnitSystem::UnitType unitType) { + using namespace Opm::unit; + + if (unitType == UnitSystem::UNIT_TYPE_FIELD) { + unitName_ = "PSIA"; + targetRateToSiConversionFactor_ = psia; // Pa -> PSI + } + else if (unitType == UnitSystem::UNIT_TYPE_METRIC) { + unitName_ = "BARSA"; + targetRateToSiConversionFactor_ = barsa; // Pa -> bar + } + else + OPM_THROW(std::logic_error, + "Unexpected unit type " << unitType); + + return unitName_; + } + + const char* unitName_; + double targetRateToSiConversionFactor_; +}; + +// no inline implementation of this since it depends on the +// WellReport type being completed first +void Summary::writeTimeStep(int writeStepIdx, + const SimulatorTimerInterface& timer, + const WellState& wellState) +{ + // create a name -> well index map + const Opm::ScheduleConstPtr schedule = eclipseState_->getSchedule(); + const auto& timeStepWells = schedule->getWells(timer.reportStepNum()); + std::map wellNameToIdxMap; + int openWellIdx = 0; + for (size_t tsWellIdx = 0; tsWellIdx < timeStepWells.size(); ++tsWellIdx) { + if (timeStepWells[tsWellIdx]->getStatus(timer.reportStepNum()) != WellCommon::SHUT ) { + wellNameToIdxMap[timeStepWells[tsWellIdx]->name()] = openWellIdx; + openWellIdx++; + } + } + + // internal view; do not move this code out of Summary! + SummaryTimeStep tstep(*this, writeStepIdx, timer); + // write all the variables + for (auto varIt = summaryReportVars_.begin(); varIt != summaryReportVars_.end(); ++varIt) { + ecl_sum_tstep_iset(tstep.ertHandle(), + smspec_node_get_params_index((*varIt)->ertHandle()), + (*varIt)->retrieveValue(writeStepIdx, timer, wellState, wellNameToIdxMap)); + } + + // write the summary file to disk + ecl_sum_fwrite(ertHandle()); +} + +void Summary::addAllWells(Opm::EclipseStateConstPtr eclipseState, + const PhaseUsage& uses) +{ + eclipseState_ = eclipseState; + auto deckUnitType = eclipseState_->getDeckUnitSystem().getType(); + + // TODO: Only create report variables that are requested with keywords + // (e.g. "WOPR") in the input files, and only for those wells that are + // mentioned in those keywords + Opm::ScheduleConstPtr schedule = eclipseState->getSchedule(); + const auto& wells = schedule->getWells(); + const int numWells = schedule->numWells(); + for (int phaseIdx = 0; phaseIdx != BlackoilPhases::MaxNumPhases; ++phaseIdx) { + const BlackoilPhases::PhaseIndex ertPhaseIdx = + static_cast (phaseIdx); + // don't bother with reporting for phases that aren't there + if (!uses.phase_used[phaseIdx]) { + continue; + } + size_t numWellTypes = sizeof(WELL_TYPES) / sizeof(WELL_TYPES[0]); + for (size_t wellTypeIdx = 0; wellTypeIdx < numWellTypes; ++wellTypeIdx) { + const WellType wellType = WELL_TYPES[wellTypeIdx]; + for (int wellIdx = 0; wellIdx != numWells; ++wellIdx) { + // W{O,G,W}{I,P}R + addWell(std::unique_ptr ( + new WellRate(*this, + eclipseState, + wells[wellIdx], + uses, + ertPhaseIdx, + wellType, + deckUnitType))); + // W{O,G,W}{I,P}T + addWell(std::unique_ptr ( + new WellTotal(*this, + eclipseState, + wells[wellIdx], + uses, + ertPhaseIdx, + wellType, + deckUnitType))); + } + } + } + + // Add BHP monitors + for (int wellIdx = 0; wellIdx != numWells; ++wellIdx) { + // In the call below: uses, phase and the well type arguments + // are not used, except to set up an index that stores the + // well indirectly. For details see the implementation of the + // WellReport constructor, and the method + // WellReport::bhp(). + BlackoilPhases::PhaseIndex ertPhaseIdx = BlackoilPhases::Liquid; + if (!uses.phase_used[BlackoilPhases::Liquid]) { + ertPhaseIdx = BlackoilPhases::Vapour; + } + addWell(std::unique_ptr ( + new WellBhp(*this, + eclipseState, + wells[wellIdx], + uses, + ertPhaseIdx, + WELL_TYPES[0], + deckUnitType))); + } +} +} // end namespace EclipseWriterDetails + + + +/** + * Convert opm-core WellType and InjectorType to eclipse welltype + */ +int EclipseWriter::eclipseWellTypeMask(WellType wellType, WellInjector::TypeEnum injectorType) +{ + int ert_well_type = IWEL_UNDOCUMENTED_ZERO; + + if (PRODUCER == wellType) { + ert_well_type = IWEL_PRODUCER; + } else if (INJECTOR == wellType) { + switch (injectorType) { + case WellInjector::WATER: + ert_well_type = IWEL_WATER_INJECTOR; + break; + case WellInjector::GAS: + ert_well_type = IWEL_GAS_INJECTOR; + break; + case WellInjector::OIL : + ert_well_type = IWEL_OIL_INJECTOR; + break; + default: + ert_well_type = IWEL_UNDOCUMENTED_ZERO; + } + } + + return ert_well_type; +} + + +/** + * Convert opm-core WellStatus to eclipse format: > 0 open, <= 0 shut + */ +int EclipseWriter::eclipseWellStatusMask(WellCommon::StatusEnum wellStatus) +{ + int well_status = 0; + + if (wellStatus == WellCommon::OPEN) { + well_status = 1; + } + return well_status; +} + + + +/** + * Convert opm-core UnitType to eclipse format: ert_ecl_unit_enum + */ +ert_ecl_unit_enum +EclipseWriter::convertUnitTypeErtEclUnitEnum(UnitSystem::UnitType unit) +{ + switch (unit) { + case UnitSystem::UNIT_TYPE_METRIC: + return ERT_ECL_METRIC_UNITS; + + case UnitSystem::UNIT_TYPE_FIELD: + return ERT_ECL_FIELD_UNITS; + + case UnitSystem::UNIT_TYPE_LAB: + return ERT_ECL_LAB_UNITS; + } + + throw std::invalid_argument("unhandled enum value"); +} + + +void EclipseWriter::writeInit(const SimulatorTimerInterface &timer) +{ + // if we don't want to write anything, this method becomes a + // no-op... + if (!enableOutput_) { + return; + } + + writeStepIdx_ = 0; + reportStepIdx_ = -1; + + EclipseWriterDetails::Init fortio(outputDir_, baseName_, /*stepIdx=*/0, eclipseState_->getIOConfigConst()); + fortio.writeHeader(numCells_, + compressedToCartesianCellIdx_, + timer, + eclipseState_, + phaseUsage_); + + IOConfigConstPtr ioConfig = eclipseState_->getIOConfigConst(); + + if (ioConfig->getWriteINITFile()) { + if (eclipseState_->hasDeckDoubleGridProperty("PERMX")) { + auto data = eclipseState_->getDoubleGridProperty("PERMX")->getData(); + EclipseWriterDetails::convertFromSiTo(data, Opm::prefix::milli * Opm::unit::darcy); + fortio.writeKeyword("PERMX", data); + } + if (eclipseState_->hasDeckDoubleGridProperty("PERMY")) { + auto data = eclipseState_->getDoubleGridProperty("PERMY")->getData(); + EclipseWriterDetails::convertFromSiTo(data, Opm::prefix::milli * Opm::unit::darcy); + fortio.writeKeyword("PERMY", data); + } + if (eclipseState_->hasDeckDoubleGridProperty("PERMZ")) { + auto data = eclipseState_->getDoubleGridProperty("PERMZ")->getData(); + EclipseWriterDetails::convertFromSiTo(data, Opm::prefix::milli * Opm::unit::darcy); + fortio.writeKeyword("PERMZ", data); + } + } + + /* Create summary object (could not do it at construction time, + since it requires knowledge of the start time). */ + { + auto eclGrid = eclipseState_->getEclipseGrid(); + auto deckUnitType = eclipseState_->getDeckUnitSystem().getType(); + bool time_in_days = true; + + if (deckUnitType == UnitSystem::UNIT_TYPE_LAB) + time_in_days = false; + + summary_.reset(new EclipseWriterDetails::Summary(outputDir_, + baseName_, + timer, + time_in_days, + eclGrid->getNX(), + eclGrid->getNY(), + eclGrid->getNZ())); + summary_->addAllWells(eclipseState_, phaseUsage_); + } +} + +// implementation of the writeTimeStep method +void EclipseWriter::writeTimeStep(const SimulatorTimerInterface& timer, + const SimulatorState& reservoirState, + const WellState& wellState, + bool isSubstep) +{ + + // if we don't want to write anything, this method becomes a + // no-op... + if (!enableOutput_) { + return; + } + + + std::vector pressure = reservoirState.pressure(); + EclipseWriterDetails::convertFromSiTo(pressure, deckToSiPressure_); + EclipseWriterDetails::restrictAndReorderToActiveCells(pressure, gridToEclipseIdx_.size(), gridToEclipseIdx_.data()); + + std::vector saturation_water; + std::vector saturation_gas; + + if (phaseUsage_.phase_used[BlackoilPhases::Aqua]) { + saturation_water = reservoirState.saturation(); + EclipseWriterDetails::extractFromStripedData(saturation_water, + /*offset=*/phaseUsage_.phase_pos[BlackoilPhases::Aqua], + /*stride=*/phaseUsage_.num_phases); + EclipseWriterDetails::restrictAndReorderToActiveCells(saturation_water, gridToEclipseIdx_.size(), gridToEclipseIdx_.data()); + } + + + if (phaseUsage_.phase_used[BlackoilPhases::Vapour]) { + saturation_gas = reservoirState.saturation(); + EclipseWriterDetails::extractFromStripedData(saturation_gas, + /*offset=*/phaseUsage_.phase_pos[BlackoilPhases::Vapour], + /*stride=*/phaseUsage_.num_phases); + EclipseWriterDetails::restrictAndReorderToActiveCells(saturation_gas, gridToEclipseIdx_.size(), gridToEclipseIdx_.data()); + } + + + + IOConfigConstPtr ioConfig = eclipseState_->getIOConfigConst(); + + + // Write restart file + if(!isSubstep && ioConfig->getWriteRestartFile(timer.reportStepNum())) + { + const size_t ncwmax = eclipseState_->getSchedule()->getMaxNumCompletionsForWells(timer.reportStepNum()); + const size_t numWells = eclipseState_->getSchedule()->numWells(timer.reportStepNum()); + std::vector wells_ptr = eclipseState_->getSchedule()->getWells(timer.reportStepNum()); + + std::vector zwell_data( numWells * Opm::EclipseWriterDetails::Restart::NZWELZ , ""); + std::vector iwell_data( numWells * Opm::EclipseWriterDetails::Restart::NIWELZ , 0 ); + std::vector icon_data( numWells * ncwmax * Opm::EclipseWriterDetails::Restart::NICONZ , 0 ); + + EclipseWriterDetails::Restart restartHandle(outputDir_, baseName_, timer.reportStepNum(), ioConfig); + + const size_t sz = wellState.bhp().size() + wellState.perfPress().size() + wellState.perfRates().size() + wellState.temperature().size() + wellState.wellRates().size(); + std::vector xwell_data( sz , 0 ); + + restartHandle.addRestartFileXwelData(wellState, xwell_data); + + for (size_t iwell = 0; iwell < wells_ptr.size(); ++iwell) { + WellConstPtr well = wells_ptr[iwell]; + { + size_t wellIwelOffset = Opm::EclipseWriterDetails::Restart::NIWELZ * iwell; + restartHandle.addRestartFileIwelData(iwell_data, timer.reportStepNum(), well , wellIwelOffset); + } + { + size_t wellIconOffset = ncwmax * Opm::EclipseWriterDetails::Restart::NICONZ * iwell; + restartHandle.addRestartFileIconData(icon_data, well->getCompletions( timer.reportStepNum() ), wellIconOffset); + } + zwell_data[ iwell * Opm::EclipseWriterDetails::Restart::NZWELZ ] = well->name().c_str(); + } + + + { + ecl_rsthead_type rsthead_data = {}; + rsthead_data.sim_time = timer.currentPosixTime(); + rsthead_data.nactive = numCells_; + rsthead_data.nx = cartesianSize_[0]; + rsthead_data.ny = cartesianSize_[1]; + rsthead_data.nz = cartesianSize_[2]; + rsthead_data.nwells = numWells; + rsthead_data.niwelz = EclipseWriterDetails::Restart::NIWELZ; + rsthead_data.nzwelz = EclipseWriterDetails::Restart::NZWELZ; + rsthead_data.niconz = EclipseWriterDetails::Restart::NICONZ; + rsthead_data.ncwmax = ncwmax; + rsthead_data.phase_sum = Opm::EclipseWriterDetails::ertPhaseMask(phaseUsage_); + rsthead_data.sim_days = Opm::unit::convert::to(timer.simulationTimeElapsed(), Opm::unit::day); //data for doubhead + + restartHandle.writeHeader(timer, + timer.reportStepNum(), + &rsthead_data); + } + + + restartHandle.add_kw(EclipseWriterDetails::Keyword(IWEL_KW, iwell_data)); + restartHandle.add_kw(EclipseWriterDetails::Keyword(ZWEL_KW, zwell_data)); + restartHandle.add_kw(EclipseWriterDetails::Keyword(OPM_XWEL, xwell_data)); + restartHandle.add_kw(EclipseWriterDetails::Keyword(ICON_KW, icon_data)); + + + EclipseWriterDetails::Solution sol(restartHandle); + sol.add(EclipseWriterDetails::Keyword("PRESSURE", pressure)); + + + // write the cell temperature + std::vector temperature = reservoirState.temperature(); + EclipseWriterDetails::convertFromSiTo(temperature, deckToSiTemperatureFactor_, deckToSiTemperatureOffset_); + EclipseWriterDetails::restrictAndReorderToActiveCells(temperature, gridToEclipseIdx_.size(), gridToEclipseIdx_.data()); + sol.add(EclipseWriterDetails::Keyword("TEMP", temperature)); + + + if (phaseUsage_.phase_used[BlackoilPhases::Aqua]) { + sol.add(EclipseWriterDetails::Keyword(EclipseWriterDetails::saturationKeywordNames[BlackoilPhases::PhaseIndex::Aqua], saturation_water)); + } + + + if (phaseUsage_.phase_used[BlackoilPhases::Vapour]) { + sol.add(EclipseWriterDetails::Keyword(EclipseWriterDetails::saturationKeywordNames[BlackoilPhases::PhaseIndex::Vapour], saturation_gas)); + } + + + const BlackoilState* blackoilState = dynamic_cast(&reservoirState); + if (blackoilState) { + // Write RS - Dissolved GOR + const std::vector& rs = blackoilState->gasoilratio(); + sol.add(EclipseWriterDetails::Keyword("RS", rs)); + + // Write RV - Volatilized oil/gas ratio + const std::vector& rv = blackoilState->rv(); + sol.add(EclipseWriterDetails::Keyword("RV", rv)); + } + } + + + //Write RFT data for current timestep to RFT file + std::shared_ptr eclipseWriteRFTHandler = std::make_shared( + compressedToCartesianCellIdx_, + numCells_, + eclipseState_->getEclipseGrid()->getCartesianSize()); + + // Write RFT file. + { + char * rft_filename = ecl_util_alloc_filename(outputDir_.c_str(), + baseName_.c_str(), + ECL_RFT_FILE, + ioConfig->getFMTOUT(), + 0); + auto unit_type = eclipseState_->getDeckUnitSystem().getType(); + ert_ecl_unit_enum ecl_unit = convertUnitTypeErtEclUnitEnum(unit_type); + std::vector wells = eclipseState_->getSchedule()->getWells(timer.reportStepNum()); + eclipseWriteRFTHandler->writeTimeStep(rft_filename, + ecl_unit, + timer, + wells, + eclipseState_->getEclipseGrid(), + pressure, + saturation_water, + saturation_gas); + free( rft_filename ); + } + + /* Summary variables (well reporting) */ + // TODO: instead of writing the header (smspec) every time, it should + // only be written when there is a change in the well configuration + // (first timestep, in practice), and reused later. but how to do this + // without keeping the complete summary in memory (which will then + // accumulate all the timesteps)? + // + // Note: The answer to the question above is still not settled, but now we do keep + // the complete summary in memory, as a member variable in the EclipseWriter class, + // instead of creating a temporary EclipseWriterDetails::Summary in this function + // every time it is called. This has been changed so that the final summary file + // will contain data from the whole simulation, instead of just the last step. + summary_->writeTimeStep(writeStepIdx_, timer, wellState); + + ++writeStepIdx_; + // store current report index + reportStepIdx_ = timer.reportStepNum(); +} + + +EclipseWriter::EclipseWriter(const parameter::ParameterGroup& params, + Opm::EclipseStateConstPtr eclipseState, + const Opm::PhaseUsage &phaseUsage, + int numCells, + const int* compressedToCartesianCellIdx) + : eclipseState_(eclipseState) + , numCells_(numCells) + , compressedToCartesianCellIdx_(compressedToCartesianCellIdx) + , gridToEclipseIdx_(numCells, int(-1) ) + , phaseUsage_(phaseUsage) +{ + const auto eclGrid = eclipseState->getEclipseGrid(); + cartesianSize_[0] = eclGrid->getNX(); + cartesianSize_[1] = eclGrid->getNY(); + cartesianSize_[2] = eclGrid->getNZ(); + + if( compressedToCartesianCellIdx ) { + // if compressedToCartesianCellIdx available then + // compute mapping to eclipse order + std::map< int , int > indexMap; + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + int cartesianCellIdx = compressedToCartesianCellIdx[cellIdx]; + indexMap[ cartesianCellIdx ] = cellIdx; + } + + int idx = 0; + for( auto it = indexMap.begin(), end = indexMap.end(); it != end; ++it ) { + gridToEclipseIdx_[ idx++ ] = (*it).second; + } + } + else { + // if not compressedToCartesianCellIdx was given use identity + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) { + gridToEclipseIdx_[ cellIdx ] = cellIdx; + } + } + + + // factor from the pressure values given in the deck to Pascals + deckToSiPressure_ = + eclipseState->getDeckUnitSystem().parse("Pressure")->getSIScaling(); + + // factor and offset from the temperature values given in the deck to Kelvin + deckToSiTemperatureFactor_ = + eclipseState->getDeckUnitSystem().parse("Temperature")->getSIScaling(); + deckToSiTemperatureOffset_ = + eclipseState->getDeckUnitSystem().parse("Temperature")->getSIOffset(); + + init(params); +} + +void EclipseWriter::init(const parameter::ParameterGroup& params) +{ + // get the base name from the name of the deck + using boost::filesystem::path; + path deckPath(params.get ("deck_filename")); + if (boost::to_upper_copy(path(deckPath.extension()).string()) == ".DATA") { + baseName_ = path(deckPath.stem()).string(); + } + else { + baseName_ = path(deckPath.filename()).string(); + } + + // make uppercase of everything (or otherwise we'll get uppercase + // of some of the files (.SMSPEC, .UNSMRY) and not others + baseName_ = boost::to_upper_copy(baseName_); + + // retrieve the value of the "output" parameter + enableOutput_ = params.getDefault("output", /*defaultValue=*/true); + + // store in current directory if not explicitly set + outputDir_ = params.getDefault("output_dir", "."); + + // set the index of the first time step written to 0... + writeStepIdx_ = 0; + reportStepIdx_ = -1; + + if (enableOutput_) { + // make sure that the output directory exists, if not try to create it + if (!boost::filesystem::exists(outputDir_)) { + std::cout << "Trying to create directory \"" << outputDir_ << "\" for the simulation output\n"; + boost::filesystem::create_directories(outputDir_); + } + + if (!boost::filesystem::is_directory(outputDir_)) { + OPM_THROW(std::runtime_error, + "The path specified as output directory '" << outputDir_ + << "' is not a directory"); + } + } +} + +// default destructor is OK, just need to be defined +EclipseWriter::~EclipseWriter() +{ } + +} // namespace Opm diff --git a/src/opm/output/eclipse/writeECLData.cpp b/src/opm/output/eclipse/writeECLData.cpp new file mode 100644 index 000000000..72129b3f2 --- /dev/null +++ b/src/opm/output/eclipse/writeECLData.cpp @@ -0,0 +1,193 @@ +/* + Copyright 2012 SINTEF ICT, Applied Mathematics. + Copyright 2012 Statoil ASA. + + 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 . +*/ + +#if HAVE_CONFIG_H +#include "config.h" +#endif + +#include +#include +#include +#include + +#include + +#ifdef HAVE_ERT // This one goes almost to the bottom of the file + +#include +#include +#include + + +namespace Opm +{ + + static ecl_kw_type * ecl_kw_wrapper( const UnstructuredGrid& grid, + 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; + } + } + + + /* + This function will write the data solution data in the DataMap + @data as an ECLIPSE restart file, in addition to the solution + fields the ECLIPSE restart file will have a minimum (hopefully + sufficient) amount of header information. + + The ECLIPSE restart files come in two varietes; unified restart + files which have all the report steps lumped together in one large + chunk and non-unified restart files which are one file for each + report step. In addition the files can be either formatted + (i.e. ASCII) or unformatted (i.e. binary). + + The writeECLData() function has two hardcoded settings: + 'file_type' and 'fmt_file' which regulate which type of files the + should be created. The extension of the files follow a convention: + + Unified, formatted : .FUNRST + Unified, unformatted : .UNRST + Multiple, formatted : .Fnnnn + Multiple, unformatted : .Xnnnn + + For the multiple files the 'nnnn' part is the report number, + formatted with '%04d' format specifier. The writeECLData() + function will use the ecl_util_alloc_filename() function to create + an ECLIPSE filename according to this conventions. + */ + + void writeECLData(const UnstructuredGrid& grid, + const DataMap& data, + const int current_step, + const double current_time, + const boost::posix_time::ptime& current_date_time, + const std::string& output_dir, + const std::string& base_name) { + + 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(output_dir.c_str() , base_name.c_str() , file_type , fmt_file , current_step ); + int phases = ECL_OIL_PHASE + ECL_WATER_PHASE; + 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 = (current_date_time - t0).total_seconds(); + + date = time_t( seconds ); + } + + if (current_step > 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_rsthead_type rsthead_data = {}; + + const int num_wells = 0; + const int niwelz = 0; + const int nzwelz = 0; + const int niconz = 0; + const int ncwmax = 0; + + rsthead_data.nx = nx; + rsthead_data.ny = ny; + rsthead_data.nz = nz; + rsthead_data.nwells = num_wells; + rsthead_data.niwelz = niwelz; + rsthead_data.nzwelz = nzwelz; + rsthead_data.niconz = niconz; + rsthead_data.ncwmax = ncwmax; + rsthead_data.nactive = nactive; + rsthead_data.phase_sum = phases; + rsthead_data.sim_time = date; + + rsthead_data.sim_days = Opm::unit::convert::to(current_time, Opm::unit::day); //Data for doubhead + + ecl_rst_file_fwrite_header( rst_file , current_step , &rsthead_data); + } + + ecl_rst_file_start_solution( rst_file ); + + { + DataMap::const_iterator i = data.find("pressure"); + if (i != data.end()) { + ecl_kw_type * pressure_kw = ecl_kw_wrapper( grid , "PRESSURE" , i->second , 0 , 1); + ecl_rst_file_add_kw( rst_file , pressure_kw ); + ecl_kw_free( pressure_kw ); + } + } + + { + DataMap::const_iterator i = data.find("saturation"); + if (i != data.end()) { + if (int(i->second->size()) != 2 * grid.number_of_cells) { + OPM_THROW(std::runtime_error, "writeECLData() requires saturation field to have two phases."); + } + ecl_kw_type * swat_kw = ecl_kw_wrapper( grid , "SWAT" , i->second , 0 , 2); + ecl_rst_file_add_kw( rst_file , swat_kw ); + ecl_kw_free( swat_kw ); + } + } + + ecl_rst_file_end_solution( rst_file ); + ecl_rst_file_close( rst_file ); + free(filename); + } +} + +#else // that is, we have not defined HAVE_ERT + +namespace Opm +{ + + void writeECLData(const UnstructuredGrid&, + const DataMap&, + const int, + const double, + const boost::posix_time::ptime&, + const std::string&, + const std::string&) + { + OPM_THROW(std::runtime_error, "Cannot call writeECLData() without ERT library support. Reconfigure opm-core with ERT support and recompile."); + } +} + +#endif diff --git a/tests/testBlackoilState3.DATA b/tests/testBlackoilState3.DATA new file mode 100644 index 000000000..27834606a --- /dev/null +++ b/tests/testBlackoilState3.DATA @@ -0,0 +1,119 @@ +RUNSPEC + +TITLE +UNIT TEST + + +START +1 NOV 1979 / + +OIL +GAS +WATER + + +METRIC + +DIMENS + 10 10 10 / + +GRID +DXV +10*10 / + +DYV +10*10 / + +DZV +10*10 / + +DEPTHZ +121*1000 / + + + + + +SCHEDULE + + +RPTSCHED + RESTART=2/ + +DATES + 1 DES 1979/ +/ + +WELSPECS + 'OP_1' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'OP_2' 'OP' 8 8 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'OP_3' 'OP' 7 7 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / +/ + +COMPDAT + 'OP_1' 9 9 1 1 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / + 'OP_1' 9 9 2 2 'OPEN' 1* 46.825 0.311 4332.346 1* 1* 'X' 22.123 / + 'OP_1' 9 9 3 3 'OPEN' 1* 51.867 0.311 4799.764 1* 1* 'X' 22.143 / + 'OP_1' 9 9 4 4 'OPEN' 1* 34.243 0.311 3169.482 1* 1* 'X' 22.166 / + 'OP_1' 9 9 5 5 'OPEN' 1* 36.435 0.311 3375.309 1* 1* 'X' 22.262 / + 'OP_1' 9 9 6 6 'OPEN' 1* 39.630 0.311 3672.067 1* 1* 'X' 22.283 / + 'OP_1' 9 9 7 7 'OPEN' 1* 33.975 0.311 3148.671 1* 1* 'X' 22.307 / + 'OP_1' 9 9 8 8 'OPEN' 1* 24.869 0.311 2305.242 1* 1* 'X' 22.329 / + 'OP_1' 9 9 9 9 'OPEN' 1* 38.301 0.311 3551.043 1* 1* 'X' 22.351 / + 'OP_1' 9 9 10 10 'OPEN' 1* 6.642 0.311 615.914 1* 1* 'X' 22.372 / + 'OP_2' 8 8 1 3 'OPEN' 1* 1.168 0.311 107.872 1* 1* 'Y' 21.925 / + 'OP_2' 8 7 3 3 'OPEN' 1* 15.071 0.311 1391.859 1* 1* 'Y' 21.920 / + 'OP_2' 8 7 3 6 'OPEN' 1* 6.242 0.311 576.458 1* 1* 'Y' 21.915 / + 'OP_3' 7 7 1 1 'OPEN' 1* 27.412 0.311 2445.337 1* 1* 'Y' 18.521 / + 'OP_3' 7 7 2 2 'OPEN' 1* 55.195 0.311 4923.842 1* 1* 'Y' 18.524 / + 'OP_3' 7 7 3 3 'OPEN' 1* 18.032 0.311 1608.615 1* 1* 'Y' 18.526 / + 'OP_3' 7 7 4 4 'OPEN' 1* 56.817 0.311 5047.177 1* 1* 'Y' 18.155 / + 'OP_3' 7 7 5 5 'OPEN' 1* 4.728 0.311 420.067 1* 1* 'Y' 18.162 / +/ + +DATES + 1 JUN 1980/ +/ + + +WELSPECS + 'OP_4' 'OP' 2 2 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'OP_5' 'OP' 5 4 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'OP_6' 'OP' 8 2 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / +/ + +COMPDAT + 'OP_4' 2 2 1 10 'OPEN' 1* 6.642 0.311 615.914 1* 1* 'X' 22.372 / + 'OP_5' 5 4 1 3 'OPEN' 1* 1.168 0.311 107.872 1* 1* 'Y' 21.925 / + 'OP_6' 8 2 1 3 'OPEN' 1* 27.412 0.311 2445.337 1* 1* 'Y' 18.521 / + 'OP_6' 8 3 3 3 'OPEN' 1* 55.195 0.311 4923.842 1* 1* 'Y' 18.524 / + 'OP_6' 8 4 3 3 'OPEN' 1* 18.032 0.311 1608.615 1* 1* 'Y' 18.526 / + 'OP_6' 8 5 3 5 'OPEN' 1* 56.817 0.311 5047.177 1* 1* 'Y' 18.155 / + 'OP_6' 8 5 5 6 'OPEN' 1* 4.728 0.311 420.067 1* 1* 'Y' 18.162 / +/ + +DATES + 1 NOV 1980/ +/ + +WELSPECS + 'WI_1' 'WI' 3 3 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'WI_2' 'WI' 3 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / + 'WI_3' 'WI' 3 6 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / +/ + +COMPDAT + 'WI_1' 3 3 1 10 'OPEN' 1* 6.642 0.311 615.914 1* 1* 'X' 22.372 / + 'WI_2' 3 9 1 7 'OPEN' 1* 1.168 0.311 107.872 1* 1* 'Y' 21.925 / + 'WI_3' 3 6 1 3 'OPEN' 1* 27.412 0.311 2445.337 1* 1* 'Y' 18.521 / +/ + + + + +DATES + 1 NOV 1982/ +/ + +END + diff --git a/tests/test_EclipseWriteRFTHandler.cpp b/tests/test_EclipseWriteRFTHandler.cpp new file mode 100755 index 000000000..80d92b89b --- /dev/null +++ b/tests/test_EclipseWriteRFTHandler.cpp @@ -0,0 +1,234 @@ +/* + Copyright 2015 Statoil ASA. + + 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" + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE EclipseRFTWriter +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include + +namespace { + +void verifyRFTFile(const std::string& rft_filename) { + + ecl_rft_file_type * new_rft_file = ecl_rft_file_alloc(rft_filename.c_str()); + std::shared_ptr rft_file; + rft_file.reset(new_rft_file, ecl_rft_file_free); + + //Get RFT node for well/time OP_1/10 OKT 2008 + time_t recording_time = util_make_datetime(0, 0, 0, 10, 10, 2008); + ecl_rft_node_type * ecl_rft_node = ecl_rft_file_get_well_time_rft(rft_file.get() , "OP_1" , recording_time); + BOOST_CHECK(ecl_rft_node_is_RFT(ecl_rft_node)); + + //Verify RFT data for completions (ijk) 9 9 1, 9 9 2 and 9 9 3 for OP_1 + const ecl_rft_cell_type * ecl_rft_cell1 = ecl_rft_node_lookup_ijk(ecl_rft_node, 8, 8, 0); + const ecl_rft_cell_type * ecl_rft_cell2 = ecl_rft_node_lookup_ijk(ecl_rft_node, 8, 8, 1); + const ecl_rft_cell_type * ecl_rft_cell3 = ecl_rft_node_lookup_ijk(ecl_rft_node, 8, 8, 2); + + BOOST_CHECK_CLOSE(ecl_rft_cell_get_pressure(ecl_rft_cell1), 210088*0.00001, 0.00001); + BOOST_CHECK_CLOSE(ecl_rft_cell_get_pressure(ecl_rft_cell2), 210188*0.00001, 0.00001); + BOOST_CHECK_CLOSE(ecl_rft_cell_get_pressure(ecl_rft_cell3), 210288*0.00001, 0.00001); + + BOOST_CHECK_EQUAL(ecl_rft_cell_get_sgas(ecl_rft_cell1), 0.0); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_sgas(ecl_rft_cell2), 0.0); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_sgas(ecl_rft_cell3), 0.0); + + BOOST_CHECK_EQUAL(ecl_rft_cell_get_swat(ecl_rft_cell1), 0.0); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_swat(ecl_rft_cell2), 0.0); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_swat(ecl_rft_cell3), 0.0); + + BOOST_CHECK_EQUAL(ecl_rft_cell_get_soil(ecl_rft_cell1), 1.0); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_soil(ecl_rft_cell2), 1.0); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_soil(ecl_rft_cell3), 1.0); + + BOOST_CHECK_EQUAL(ecl_rft_cell_get_depth(ecl_rft_cell1), (0.250 + (0.250/2))); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_depth(ecl_rft_cell2), (2*0.250 + (0.250/2))); + BOOST_CHECK_EQUAL(ecl_rft_cell_get_depth(ecl_rft_cell3), (3*0.250 + (0.250/2))); +} + + + + +Opm::DeckConstPtr createDeck(const std::string& input_str) { + Opm::ParserPtr parser = std::make_shared(); + Opm::DeckConstPtr deck = parser->parseString(input_str , Opm::ParseMode()); + return deck; +} + + +std::shared_ptr createWellState(std::shared_ptr blackoilState) +{ + std::shared_ptr wellState = std::make_shared(); + wellState->init(0, *blackoilState); + return wellState; +} + + + +std::shared_ptr createBlackoilState(int timeStepIdx, std::shared_ptr ourFineGridManagerPtr) +{ + const UnstructuredGrid &ourFinerUnstructuredGrid = *ourFineGridManagerPtr->c_grid(); + + std::shared_ptr blackoilState = std::make_shared(); + blackoilState->init(ourFinerUnstructuredGrid, 3); + + size_t numCells = ourFinerUnstructuredGrid.number_of_cells; + + auto &pressure = blackoilState->pressure(); + for (size_t cellIdx = 0; cellIdx < numCells; ++cellIdx) { + pressure[cellIdx] = timeStepIdx*1e5 + 1e4 + cellIdx; + } + return blackoilState; +} + + + +std::shared_ptr createEclipseWriter(std::shared_ptr deck, + std::shared_ptr eclipseState, + std::shared_ptr ourFineGridManagerPtr, + const int * compressedToCartesianCellIdx) +{ + Opm::parameter::ParameterGroup params; + params.insertParameter("deck_filename", "testcase.data"); + + Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck); + + const UnstructuredGrid &ourFinerUnstructuredGrid = *ourFineGridManagerPtr->c_grid(); + + std::shared_ptr eclipseWriter = std::make_shared(params, + eclipseState, + phaseUsage, + ourFinerUnstructuredGrid.number_of_cells, + compressedToCartesianCellIdx); + + return eclipseWriter; +} + +} + +BOOST_AUTO_TEST_CASE(test_EclipseWriterRFTHandler) +{ + const std::string& deckString = + "RUNSPEC\n" + "OIL\n" + "GAS\n" + "WATER\n" + "DIMENS\n" + " 10 10 10 /\n" + "GRID\n" + "DXV\n" + "10*0.25 /\n" + "DYV\n" + "10*0.25 /\n" + "DZV\n" + "10*0.25 /\n" + "TOPS\n" + "100*0.25 /\n" + "\n" + "START -- 0 \n" + "1 NOV 1979 / \n" + "SCHEDULE\n" + "DATES -- 1\n" + " 1 DES 1979/ \n" + "/\n" + "WELSPECS\n" + " 'OP_1' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + " 'OP_2' 'OP' 4 4 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + "/\n" + "COMPDAT\n" + " 'OP_1' 9 9 1 1 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + " 'OP_1' 9 9 2 2 'OPEN' 1* 46.825 0.311 4332.346 1* 1* 'X' 22.123 / \n" + " 'OP_1' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + " 'OP_2' 4 4 4 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + "DATES -- 2\n" + " 10 OKT 2008 / \n" + "/\n" + "WRFT \n" + "/ \n" + "WELOPEN\n" + " 'OP_1' OPEN / \n" + " 'OP_2' OPEN / \n" + "/\n" + "DATES -- 3\n" + " 10 NOV 2008 / \n" + "/\n"; + + + + test_work_area_type * new_ptr = test_work_area_alloc("test_EclipseWriterRFTHandler"); + std::shared_ptr test_area; + test_area.reset(new_ptr, test_work_area_free); + + std::shared_ptr deck = createDeck(deckString); + std::shared_ptr eclipseState = std::make_shared(deck , Opm::ParseMode()); + + std::shared_ptr simulatorTimer = std::make_shared(); + simulatorTimer->init(eclipseState->getSchedule()->getTimeMap()); + + std::shared_ptr ourFineGridManagerPtr = std::make_shared(eclipseState->getEclipseGrid()); + const UnstructuredGrid &ourFinerUnstructuredGrid = *ourFineGridManagerPtr->c_grid(); + const int* compressedToCartesianCellIdx = Opm::UgGridHelpers::globalCell(ourFinerUnstructuredGrid); + + std::shared_ptr eclipseWriter = createEclipseWriter(deck, + eclipseState, + ourFineGridManagerPtr, + compressedToCartesianCellIdx); + eclipseWriter->writeInit(*simulatorTimer); + + + for (; simulatorTimer->currentStepNum() < simulatorTimer->numSteps(); ++ (*simulatorTimer)) { + std::shared_ptr blackoilState2 = createBlackoilState(simulatorTimer->currentStepNum(),ourFineGridManagerPtr); + std::shared_ptr wellState = createWellState(blackoilState2); + eclipseWriter->writeTimeStep(*simulatorTimer, *blackoilState2, *wellState, false); + } + + std::string cwd(test_work_area_get_cwd(test_area.get())); + std::string rft_filename = cwd + "/TESTCASE.RFT"; + verifyRFTFile(rft_filename); + +} + + + diff --git a/tests/test_EclipseWriter.cpp b/tests/test_EclipseWriter.cpp new file mode 100644 index 000000000..020ad18f6 --- /dev/null +++ b/tests/test_EclipseWriter.cpp @@ -0,0 +1,421 @@ +/* + Copyright 2014 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 . +*/ +#include "config.h" + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE EclipseWriter +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +// ERT stuff +#include +#include +#include + +#include + +std::shared_ptr eclWriter; +std::shared_ptr simTimer; +std::shared_ptr deck; +std::shared_ptr eclipseState; +std::shared_ptr ourFineGridManagerPtr; +std::shared_ptr blackoilState; +std::shared_ptr wellState; + +void createEclipseWriter(const char *deckString) +{ + Opm::ParseMode parseMode; + Opm::ParserConstPtr parser(new Opm::Parser()); + deck = parser->parseString(deckString, parseMode); + + Opm::parameter::ParameterGroup params; + params.insertParameter("deck_filename", "foo.data"); + + eclipseState.reset(new Opm::EclipseState(deck , parseMode)); + + auto eclGrid = eclipseState->getEclipseGrid(); + BOOST_CHECK(eclGrid->getNX() == 3); + BOOST_CHECK(eclGrid->getNY() == 3); + BOOST_CHECK(eclGrid->getNZ() == 3); + BOOST_CHECK(eclGrid->getCartesianSize() == 3*3*3); + + simTimer.reset(new Opm::SimulatorTimer()); + simTimer->init(eclipseState->getSchedule()->getTimeMap()); + + // also create an UnstructuredGrid (required to create a BlackoilState) + Opm::EclipseGridConstPtr constEclGrid(eclGrid); + ourFineGridManagerPtr.reset(new Opm::GridManager(constEclGrid)); + + const UnstructuredGrid &ourFinerUnstructuredGrid = *ourFineGridManagerPtr->c_grid(); + BOOST_CHECK(ourFinerUnstructuredGrid.cartdims[0] == 3); + BOOST_CHECK(ourFinerUnstructuredGrid.cartdims[1] == 3); + BOOST_CHECK(ourFinerUnstructuredGrid.cartdims[2] == 3); + + BOOST_CHECK(ourFinerUnstructuredGrid.number_of_cells == 3*3*3); + + Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck); + eclWriter.reset(new Opm::EclipseWriter(params, + eclipseState, + phaseUsage, + ourFinerUnstructuredGrid.number_of_cells, + 0)); + + // this check is disabled so far, because UnstructuredGrid uses some weird definition + // of the term "face". For this grid, "number_of_faces" is 108 which is + // 2*6*numCells... + //BOOST_CHECK(ourFinerUnstructuredGrid.number_of_faces == 4*4*4); + + int numCells = ourFinerUnstructuredGrid.number_of_cells; + for (int cellIdx = 0; cellIdx < numCells; ++cellIdx) + BOOST_CHECK(ourFinerUnstructuredGrid.global_cell[cellIdx] == cellIdx); +} + +void createBlackoilState(int timeStepIdx) +{ + // allocate a new BlackoilState object + const UnstructuredGrid &ourFinerUnstructuredGrid = *ourFineGridManagerPtr->c_grid(); + blackoilState.reset(new Opm::BlackoilState); + blackoilState->init(ourFinerUnstructuredGrid, 3); + + size_t numCells = ourFinerUnstructuredGrid.number_of_cells; + size_t numFaces = ourFinerUnstructuredGrid.number_of_faces; + + BOOST_CHECK(blackoilState->pressure().size() == numCells); + BOOST_CHECK(blackoilState->facepressure().size() == numFaces); + BOOST_CHECK(blackoilState->faceflux().size() == numFaces); + BOOST_CHECK(blackoilState->saturation().size() == numCells*3); + BOOST_CHECK(blackoilState->gasoilratio().size() == numCells); + BOOST_CHECK(blackoilState->rv().size() == numCells); + + // this check is disabled because BlackoilState does not seem to allocate memory for + // this field. This means that it is probably unused and unneeded. + //BOOST_CHECK(blackoilState->surfacevol().size() == numCells*3); + + // fill the state object with some data. The fun with this class is that it does not + // exhibit a proper c++ way to do this (i.e., getter + setter methods). Instead + // references to the arrays must be retrieved from the object and manipulated + // directly. Don't try to call resize() or anything else which is not politically + // correct on them! + auto &pressure = blackoilState->pressure(); + auto &facepressure = blackoilState->facepressure(); + auto &faceflux = blackoilState->faceflux(); + auto &saturation = blackoilState->saturation(); + auto &gasoilratio = blackoilState->gasoilratio(); + auto &rv = blackoilState->rv(); + for (size_t cellIdx = 0; cellIdx < numCells; ++cellIdx) { + pressure[cellIdx] = timeStepIdx*1e5 + 1e4 + cellIdx; + + // set the phase saturations. Some fun with direct index manipulation is to be + // had... + saturation[3*cellIdx + 0] = timeStepIdx*1e5 +2.1e4 + cellIdx; // oil + saturation[3*cellIdx + 1] = timeStepIdx*1e5 +2.2e4 + cellIdx; // gas + saturation[3*cellIdx + 2] = timeStepIdx*1e5 +2.3e4 + cellIdx; // water + + // oil vaporization factor + rv[cellIdx] = timeStepIdx*1e5 +3e4 + cellIdx; + + // gas dissolution factor + gasoilratio[cellIdx] = timeStepIdx*1e5 + 4e4 + cellIdx; + } + + // face specific data + for (size_t faceIdx = 0; faceIdx < numFaces; ++faceIdx) { + facepressure[faceIdx] = timeStepIdx*1e5 + 5e4 + faceIdx; + faceflux[faceIdx] = timeStepIdx*1e5 + 6e4 + faceIdx; + } +} + +void createWellState(int /*timeStepIdx*/) +{ + // allocate a new BlackoilState object + wellState.reset(new Opm::WellState); + wellState->init(0, *blackoilState); +} + +void getErtData(ecl_kw_type *eclKeyword, std::vector &data) +{ + size_t kwSize = ecl_kw_get_size(eclKeyword); + float* ertData = static_cast(ecl_kw_iget_ptr(eclKeyword, 0)); + + data.resize(kwSize); + std::copy(ertData, ertData + kwSize, data.begin()); +} + +void getErtData(ecl_kw_type *eclKeyword, std::vector &data) +{ + size_t kwSize = ecl_kw_get_size(eclKeyword); + int* ertData = static_cast(ecl_kw_iget_ptr(eclKeyword, 0)); + + data.resize(kwSize); + std::copy(ertData, ertData + kwSize, data.begin()); +} + +void compareErtData(const std::vector &src, const std::vector &dst, double tolerance) +{ + BOOST_CHECK_EQUAL(src.size(), dst.size()); + if (src.size() != dst.size()) + return; + + for (size_t i = 0; i < src.size(); ++i) + BOOST_CHECK_CLOSE(src[i], dst[i], tolerance); +} + +void compareErtData(const std::vector &src, const std::vector &dst) +{ + BOOST_CHECK_EQUAL(src.size(), dst.size()); + if (src.size() != dst.size()) + return; + + for (size_t i = 0; i < src.size(); ++i) + BOOST_CHECK_EQUAL(src[i], dst[i]); +} + +void checkEgridFile() +{ + size_t numCells = ourFineGridManagerPtr->c_grid()->number_of_cells; + + // use ERT directly to inspect the EGRID file produced by EclipseWriter + auto egridFile = fortio_open_reader("FOO.EGRID", /*isFormated=*/0, ECL_ENDIAN_FLIP); + + auto eclGrid = eclipseState->getEclipseGrid(); + + ecl_kw_type *eclKeyword; + // yes, that's an assignment! + while ((eclKeyword = ecl_kw_fread_alloc(egridFile))) { + std::string keywordName(ecl_kw_get_header(eclKeyword)); + if (keywordName == "COORD") { + std::vector sourceData, resultData; + eclGrid->exportCOORD(sourceData); + getErtData(eclKeyword, resultData); + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-6); + } + else if (keywordName == "ZCORN") { + std::vector sourceData, resultData; + eclGrid->exportZCORN(sourceData); + getErtData(eclKeyword, resultData); + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-6); + } + else if (keywordName == "ACTNUM") { + std::vector sourceData, resultData; + eclGrid->exportACTNUM(sourceData); + getErtData(eclKeyword, resultData); + + if (resultData.size() == numCells && sourceData.size() == 0) { + sourceData.resize(numCells); + std::fill(sourceData.begin(), sourceData.end(), 1); + } + + compareErtData(sourceData, resultData); + } + + ecl_kw_free(eclKeyword); + } + + fortio_fclose(egridFile); +} + +void checkInitFile() +{ + // use ERT directly to inspect the INIT file produced by EclipseWriter + auto initFile = fortio_open_reader("FOO.INIT", /*isFormated=*/0, ECL_ENDIAN_FLIP); + + ecl_kw_type *eclKeyword; + // yes, that's an assignment! + while ((eclKeyword = ecl_kw_fread_alloc(initFile))) { + std::string keywordName(ecl_kw_get_header(eclKeyword)); + + if (keywordName == "PORO") { + const std::vector &sourceData = deck->getKeyword("PORO").getSIDoubleData(); + std::vector resultData; + getErtData(eclKeyword, resultData); + + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-4); + } + + if (keywordName == "PERMX") { + std::vector sourceData = deck->getKeyword("PERMX").getSIDoubleData(); + std::vector resultData; + getErtData(eclKeyword, resultData); + + // convert the data from ERT from Field to SI units (mD to m^2) + for (size_t i = 0; i < resultData.size(); ++i) { + resultData[i] *= 9.869233e-16; + } + + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-4); + } + + ecl_kw_free(eclKeyword); + } + + fortio_fclose(initFile); +} + +void checkRestartFile(int timeStepIdx) +{ + size_t numCells = ourFineGridManagerPtr->c_grid()->number_of_cells; + + Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck); + int numActivePhases = phaseUsage.num_phases; + int waterPhaseIdx = phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua]; + int gasPhaseIdx = phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour]; + + for (int i = 0; i <= timeStepIdx; ++i) { + createBlackoilState(i); + + // use ERT directly to inspect the restart file produced by EclipseWriter + auto rstFile = fortio_open_reader("FOO.UNRST", /*isFormated=*/0, ECL_ENDIAN_FLIP); + + int curSeqnum = -1; + ecl_kw_type *eclKeyword; + // yes, that's an assignment! + while ((eclKeyword = ecl_kw_fread_alloc(rstFile))) { + std::string keywordName(ecl_kw_get_header(eclKeyword)); + + if (keywordName == "SEQNUM") { + curSeqnum = *static_cast(ecl_kw_iget_ptr(eclKeyword, 0)); + } + if (curSeqnum != i) + continue; + + if (keywordName == "PRESSURE") { + std::vector sourceData = blackoilState->pressure(); + std::vector resultData; + getErtData(eclKeyword, resultData); + + // convert the data from ERT from Metric to SI units (bar to Pa) + for (size_t ii = 0; ii < resultData.size(); ++ii) { + resultData[ii] *= 1e5; + } + + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-4); + } + + if (keywordName == "SWAT") { + std::vector sourceData; + std::vector resultData; + getErtData(eclKeyword, resultData); + + // extract the water saturation from the black-oil state + sourceData.resize(numCells); + for (size_t ii = 0; ii < sourceData.size(); ++ii) { + // again, fun with direct index manipulation... + sourceData[ii] = blackoilState->saturation()[ii*numActivePhases + waterPhaseIdx]; + } + + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-4); + } + + if (keywordName == "SGAS") { + std::vector sourceData; + std::vector resultData; + getErtData(eclKeyword, resultData); + + // extract the water saturation from the black-oil state + sourceData.resize(numCells); + for (size_t ii = 0; ii < sourceData.size(); ++ii) { + // again, fun with direct index manipulation... + sourceData[ii] = blackoilState->saturation()[ii*numActivePhases + gasPhaseIdx]; + } + + compareErtData(sourceData, resultData, /*percentTolerance=*/1e-4); + } + } + + fortio_fclose(rstFile); + } +} + +void checkSummaryFile(int /*timeStepIdx*/) +{ + // TODO +} + +BOOST_AUTO_TEST_CASE(EclipseWriterIntegration) +{ + const char *deckString = + "RUNSPEC\n" + "INIT\n" + "UNIFOUT\n" + "OIL\n" + "GAS\n" + "WATER\n" + "METRIC\n" + "DIMENS\n" + "3 3 3/\n" + "GRID\n" + "DXV\n" + "1.0 2.0 3.0 /\n" + "DYV\n" + "4.0 5.0 6.0 /\n" + "DZV\n" + "7.0 8.0 9.0 /\n" + "TOPS\n" + "9*100 /\n" + "PROPS\n" + "PORO\n" + "27*0.3 /\n" + "PERMX\n" + "27*1 /\n" + "SOLUTION\n" + "RPTRST\n" + "BASIC=2\n" + "/\n" + "SCHEDULE\n" + "TSTEP\n" + "1.0 2.0 3.0 4.0 /\n" + "WELSPECS\n" + "'INJ' 'G' 1 1 2000 'GAS' /\n" + "'PROD' 'G' 3 3 1000 'OIL' /\n" + "/\n"; + + createEclipseWriter(deckString); + + eclWriter->writeInit(*simTimer); + + checkEgridFile(); + checkInitFile(); + + for (; simTimer->currentStepNum() < simTimer->numSteps(); ++ (*simTimer)) { + createBlackoilState(simTimer->currentStepNum()); + createWellState(simTimer->currentStepNum()); + eclWriter->writeTimeStep(*simTimer, *blackoilState, *wellState, false); + checkRestartFile(simTimer->currentStepNum()); + checkSummaryFile(simTimer->currentStepNum()); + } +} diff --git a/tests/test_writeReadRestartFile.cpp b/tests/test_writeReadRestartFile.cpp new file mode 100644 index 000000000..61584cc70 --- /dev/null +++ b/tests/test_writeReadRestartFile.cpp @@ -0,0 +1,348 @@ + +/* + Copyright 2014 Statoil IT + 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" + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE EclipseWriter +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// ERT stuff +#include +#include +#include +#include +#include +#include + +#include + + +std::string input = + "RUNSPEC\n" + "OIL\n" + "GAS\n" + "WATER\n" + "DISGAS\n" + "VAPOIL\n" + "UNIFOUT\n" + "UNIFIN\n" + "DIMENS\n" + " 10 10 10 /\n" + + "GRID\n" + "DXV\n" + "10*0.25 /\n" + "DYV\n" + "10*0.25 /\n" + "DZV\n" + "10*0.25 /\n" + "TOPS\n" + "100*0.25 /\n" + "\n" + + "SOLUTION\n" + "RESTART\n" + "TESTWELLSTATE 1/\n" + "\n" + + "START -- 0 \n" + "1 NOV 1979 / \n" + + "SCHEDULE\n" + "SKIPREST\n" + "RPTRST\n" + "BASIC=1\n" + "/\n" + "DATES -- 1\n" + " 10 OKT 2008 / \n" + "/\n" + "WELSPECS\n" + " 'OP_1' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + " 'OP_2' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + "/\n" + "COMPDAT\n" + " 'OP_1' 9 9 1 1 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + " 'OP_2' 9 9 2 2 'OPEN' 1* 46.825 0.311 4332.346 1* 1* 'X' 22.123 / \n" + " 'OP_1' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + "WCONPROD\n" + "'OP_1' 'OPEN' 'ORAT' 20000 4* 1000 /\n" + "/\n" + "WCONINJE\n" + "'OP_2' 'GAS' 'OPEN' 'RATE' 100 200 400 /\n" + "/\n" + + "DATES -- 2\n" + " 20 JAN 2011 / \n" + "/\n" + "WELSPECS\n" + " 'OP_3' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + "/\n" + "COMPDAT\n" + " 'OP_3' 9 9 1 1 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + "WCONPROD\n" + "'OP_3' 'OPEN' 'ORAT' 20000 4* 1000 /\n" + "/\n" + + "DATES -- 3\n" + " 15 JUN 2013 / \n" + "/\n" + "COMPDAT\n" + " 'OP_2' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + " 'OP_1' 9 9 7 7 'SHUT' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + + "DATES -- 4\n" + " 22 APR 2014 / \n" + "/\n" + "WELSPECS\n" + " 'OP_4' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + "/\n" + "COMPDAT\n" + " 'OP_4' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + " 'OP_3' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + "WCONPROD\n" + "'OP_4' 'OPEN' 'ORAT' 20000 4* 1000 /\n" + "/\n" + + "DATES -- 5\n" + " 30 AUG 2014 / \n" + "/\n" + "WELSPECS\n" + " 'OP_5' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + "/\n" + "COMPDAT\n" + " 'OP_5' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + "WCONPROD\n" + "'OP_5' 'OPEN' 'ORAT' 20000 4* 1000 /\n" + "/\n" + + "DATES -- 6\n" + " 15 SEP 2014 / \n" + "/\n" + "WCONPROD\n" + "'OP_3' 'SHUT' 'ORAT' 20000 4* 1000 /\n" + "/\n" + + "DATES -- 7\n" + " 9 OCT 2014 / \n" + "/\n" + "WELSPECS\n" + " 'OP_6' 'OP' 9 9 1* 'OIL' 1* 1* 1* 1* 1* 1* 1* / \n" + "/\n" + "COMPDAT\n" + " 'OP_6' 9 9 3 9 'OPEN' 1* 32.948 0.311 3047.839 1* 1* 'X' 22.100 / \n" + "/\n" + "WCONPROD\n" + "'OP_6' 'OPEN' 'ORAT' 20000 4* 1000 /\n" + "/\n" + "TSTEP -- 8\n" + "10 /" + "/\n"; + +std::shared_ptr createBlackOilState(Opm::EclipseGridConstPtr eclGrid) { + + std::shared_ptr ourFineGridManagerPtr(new Opm::GridManager(eclGrid)); + std::shared_ptr blackoilState(new Opm::BlackoilState); + blackoilState->init(*ourFineGridManagerPtr->c_grid(), 3); + + return blackoilState; +} + +Opm::EclipseWriterPtr createEclipseWriter(Opm::DeckConstPtr deck, + Opm::EclipseStatePtr eclipseState, + std::string& eclipse_data_filename) { + + Opm::parameter::ParameterGroup params; + params.insertParameter("deck_filename", eclipse_data_filename); + + const Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck); + + Opm::EclipseWriterPtr eclWriter(new Opm::EclipseWriter(params, + eclipseState, + phaseUsage, + eclipseState->getEclipseGrid()->getCartesianSize(), + 0)); + return eclWriter; +} + +void setValuesInWellState(std::shared_ptr wellState){ + wellState->bhp()[0] = 1.23; + wellState->bhp()[1] = 2.34; + + wellState->temperature()[0] = 3.45; + wellState->temperature()[1] = 4.56; + + wellState->wellRates()[0] = 5.67; + wellState->wellRates()[1] = 6.78; + wellState->wellRates()[2] = 7.89; + wellState->wellRates()[3] = 8.90; + wellState->wellRates()[4] = 9.01; + wellState->wellRates()[5] = 10.12; + + wellState->perfPress()[0] = 20.41; + wellState->perfPress()[1] = 21.19; + wellState->perfPress()[2] = 22.41; + wellState->perfPress()[3] = 23.19; + wellState->perfPress()[4] = 24.41; + wellState->perfPress()[5] = 25.19; + wellState->perfPress()[6] = 26.41; + wellState->perfPress()[7] = 27.19; + wellState->perfPress()[8] = 28.41; + + wellState->perfRates()[0] = 30.45; + wellState->perfRates()[1] = 31.19; + wellState->perfRates()[2] = 32.45; + wellState->perfRates()[3] = 33.19; + wellState->perfRates()[4] = 34.45; + wellState->perfRates()[5] = 35.19; + wellState->perfRates()[6] = 36.45; + wellState->perfRates()[7] = 37.19; + wellState->perfRates()[8] = 38.45; +} + +BOOST_AUTO_TEST_CASE(EclipseReadWriteWellStateData) +{ + std::string eclipse_data_filename = "TestWellState.DATA"; + test_work_area_type * test_area = test_work_area_alloc("EclipseReadWriteWellStateData"); + + Opm::Parser parser; + Opm::ParseMode parseMode; + Opm::DeckConstPtr deck = parser.parseString(input, parseMode); + Opm::EclipseStatePtr eclipseState(new Opm::EclipseState(deck , parseMode)); + Opm::EclipseWriterPtr eclipseWriter = createEclipseWriter(deck, eclipseState, eclipse_data_filename); + + std::shared_ptr simTimer( new Opm::SimulatorTimer() ); + simTimer->init(eclipseState->getSchedule()->getTimeMap()); + eclipseWriter->writeInit(*simTimer); + std::shared_ptr wellState(new Opm::WellState()); + Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck); + + Opm::GridManager gridManager(deck); + Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL); + const Wells* wells = wellsManager.c_wells(); + std::shared_ptr blackoilState = createBlackOilState(eclipseState->getEclipseGrid()); + wellState->init(wells, *blackoilState); + + //Set test data for pressure + std::vector& pressure = blackoilState->pressure(); + for (std::vector::iterator iter = pressure.begin(); iter != pressure.end(); ++iter) { + *iter = 6.0; + } + + //Set test data for temperature + std::vector& temperature = blackoilState->temperature(); + for (std::vector::iterator iter = temperature.begin(); iter != temperature.end(); ++iter) { + *iter = 7.0; + } + + //Set test data for saturation water + std::vector swatdata(1000, 8); + Opm::EclipseIOUtil::addToStripedData(swatdata, blackoilState->saturation(), phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua], phaseUsage.num_phases); + + //Set test data for saturation gas + std::vector sgasdata(1000, 9); + Opm::EclipseIOUtil::addToStripedData(sgasdata, blackoilState->saturation(), phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour], phaseUsage.num_phases); + + // Set test data for rs + double rs = 300.0; + std::vector& rs_vec = blackoilState->gasoilratio(); + for (std::vector::iterator rs_iter = rs_vec.begin(); rs_iter != rs_vec.end(); ++ rs_iter) { + *rs_iter = rs; + rs = rs + 1.0; + } + + // Set testdata for rv + double rv = 400.0; + std::vector& rv_vec = blackoilState->rv(); + for (std::vector::iterator rv_iter = rv_vec.begin(); rv_iter != rv_vec.end(); ++rv_iter) { + *rv_iter = rv; + rv = rv + 1.0; + } + + setValuesInWellState(wellState); + simTimer->setCurrentStepNum(1); + eclipseWriter->writeTimeStep(*simTimer, *blackoilState, *wellState , false); + + std::shared_ptr wellStateRestored(new Opm::WellState()); + wellStateRestored->init(wells, *blackoilState); + + //Read and verify OPM XWEL data, and solution data: pressure, temperature, saturation data, rs and rv + std::shared_ptr blackoilStateRestored = createBlackOilState(eclipseState->getEclipseGrid()); + Opm::init_from_restart_file(eclipseState, Opm::UgGridHelpers::numCells(*gridManager.c_grid()), phaseUsage, *blackoilStateRestored, *wellStateRestored); + + BOOST_CHECK_EQUAL_COLLECTIONS(wellState->bhp().begin(), wellState->bhp().end(), wellStateRestored->bhp().begin(), wellStateRestored->bhp().end()); + BOOST_CHECK_EQUAL_COLLECTIONS(wellState->temperature().begin(), wellState->temperature().end(), wellStateRestored->temperature().begin(), wellStateRestored->temperature().end()); + BOOST_CHECK_EQUAL_COLLECTIONS(wellState->wellRates().begin(), wellState->wellRates().end(), wellStateRestored->wellRates().begin(), wellStateRestored->wellRates().end()); + BOOST_CHECK_EQUAL_COLLECTIONS(wellState->perfRates().begin(), wellState->perfRates().end(), wellStateRestored->perfRates().begin(), wellStateRestored->perfRates().end()); + BOOST_CHECK_EQUAL_COLLECTIONS(wellState->perfPress().begin(), wellState->perfPress().end(), wellStateRestored->perfPress().begin(), wellStateRestored->perfPress().end()); + + std::vector swat_restored; + std::vector swat; + std::vector sgas_restored; + std::vector sgas; + Opm::EclipseIOUtil::extractFromStripedData(blackoilStateRestored->saturation(), swat_restored, phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua], phaseUsage.num_phases); + Opm::EclipseIOUtil::extractFromStripedData(blackoilState->saturation(), swat, phaseUsage.phase_pos[Opm::BlackoilPhases::Aqua], phaseUsage.num_phases); + Opm::EclipseIOUtil::extractFromStripedData(blackoilStateRestored->saturation(), sgas_restored, phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour], phaseUsage.num_phases); + Opm::EclipseIOUtil::extractFromStripedData(blackoilState->saturation(), sgas, phaseUsage.phase_pos[Opm::BlackoilPhases::Vapour], phaseUsage.num_phases); + + for (size_t cellindex = 0; cellindex < 1000; ++cellindex) { + BOOST_CHECK_CLOSE(blackoilState->pressure()[cellindex], blackoilStateRestored->pressure()[cellindex], 0.00001); + BOOST_CHECK_CLOSE(blackoilState->temperature()[cellindex], blackoilStateRestored->temperature()[cellindex], 0.00001); + BOOST_CHECK_CLOSE(swat[cellindex], swat_restored[cellindex], 0.00001); + BOOST_CHECK_CLOSE(sgas[cellindex], sgas_restored[cellindex], 0.00001); + } + + + for (size_t cellindex = 0; cellindex < 1000; ++cellindex) { + BOOST_CHECK_CLOSE(blackoilState->gasoilratio()[cellindex], blackoilStateRestored->gasoilratio()[cellindex], 0.0000001); + BOOST_CHECK_CLOSE(blackoilState->rv()[cellindex], blackoilStateRestored->rv()[cellindex], 0.0000001); + } + + test_work_area_free(test_area); +} diff --git a/tests/test_writenumwells.cpp b/tests/test_writenumwells.cpp new file mode 100644 index 000000000..0d22be742 --- /dev/null +++ b/tests/test_writenumwells.cpp @@ -0,0 +1,197 @@ + +/* + Copyright 2014 Statoil IT + 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" + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define BOOST_TEST_MODULE EclipseWriter +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include + +// ERT stuff +#include +#include +#include +#include +#include +#include + +#include + + +void verifyWellState(const std::string& rst_filename, + Opm::EclipseGridConstPtr ecl_grid, + Opm::ScheduleConstPtr schedule) { + + well_info_type * well_info = well_info_alloc(ecl_grid->c_ptr()); + well_info_load_rstfile(well_info, rst_filename.c_str(), false); + + //Verify numwells + int numwells = well_info_get_num_wells(well_info); + BOOST_CHECK(numwells == (int)schedule->numWells()); + + std::vector wells = schedule->getWells(); + + for (int i = 0; i < numwells; ++i) { + + //Verify wellnames + const char * wellname = well_info_iget_well_name(well_info, i); + Opm::WellConstPtr well = wells.at(i); + BOOST_CHECK(wellname == well->name()); + + // Verify well-head position data + well_ts_type * well_ts = well_info_get_ts(well_info , wellname); + well_state_type * well_state = well_ts_iget_state(well_ts, 0); + const well_conn_type * well_head = well_state_get_wellhead(well_state, ECL_GRID_GLOBAL_GRID); + BOOST_CHECK(well_conn_get_i(well_head) == well->getHeadI()); + BOOST_CHECK(well_conn_get_j(well_head) == well->getHeadJ()); + + for (int j = 0; j < well_ts_get_size(well_ts); ++j) { + well_state = well_ts_iget_state(well_ts, j); + + //Verify welltype + int ert_well_type = well_state_get_type(well_state); + WellType welltype = well->isProducer(j) ? PRODUCER : INJECTOR; + Opm::WellInjector::TypeEnum injectortype = well->getInjectionProperties(j).injectorType; + int ecl_converted_welltype = Opm::EclipseWriter::eclipseWellTypeMask(welltype, injectortype); + int ert_converted_welltype = well_state_translate_ecl_type_int(ecl_converted_welltype); + BOOST_CHECK(ert_well_type == ert_converted_welltype); + + //Verify wellstatus + int ert_well_status = well_state_is_open(well_state) ? 1 : 0; + + Opm::WellCommon::StatusEnum status = well->getStatus(j); + int wellstatus = Opm::EclipseWriter::eclipseWellStatusMask(status); + + BOOST_CHECK(ert_well_status == wellstatus); + + //Verify number of completion connections + const well_conn_collection_type * well_connections = well_state_get_global_connections( well_state ); + size_t num_wellconnections = well_conn_collection_get_size(well_connections); + + int report_nr = well_state_get_report_nr(well_state); + Opm::CompletionSetConstPtr completions_set = well->getCompletions((size_t)report_nr); + + BOOST_CHECK(num_wellconnections == completions_set->size()); + + //Verify coordinates for each completion connection + for (size_t k = 0; k < num_wellconnections; ++k) { + const well_conn_type * well_connection = well_conn_collection_iget_const(well_connections , k); + + Opm::CompletionConstPtr completion = completions_set->get(k); + + BOOST_CHECK(well_conn_get_i(well_connection) == completion->getI()); + BOOST_CHECK(well_conn_get_j(well_connection) == completion->getJ()); + BOOST_CHECK(well_conn_get_k(well_connection) == completion->getK()); + } + } + } + + well_info_free(well_info); +} + + +std::shared_ptr createBlackOilState(Opm::EclipseGridConstPtr eclGrid) { + + std::shared_ptr ourFineGridManagerPtr(new Opm::GridManager(eclGrid)); + std::shared_ptr blackoilState(new Opm::BlackoilState); + blackoilState->init(*ourFineGridManagerPtr->c_grid(), 3); + + return blackoilState; +} + + +Opm::DeckConstPtr createDeck(const std::string& eclipse_data_filename) { + Opm::ParserPtr parser(new Opm::Parser()); + Opm::DeckConstPtr deck = parser->parseFile(eclipse_data_filename , Opm::ParseMode()); + + return deck; +} + + +Opm::EclipseWriterPtr createEclipseWriter(Opm::DeckConstPtr deck, + Opm::EclipseStatePtr eclipseState, + std::string& eclipse_data_filename) { + + Opm::parameter::ParameterGroup params; + params.insertParameter("deck_filename", eclipse_data_filename); + + const Opm::PhaseUsage phaseUsage = Opm::phaseUsageFromDeck(deck); + + Opm::EclipseWriterPtr eclWriter(new Opm::EclipseWriter(params, + eclipseState, + phaseUsage, + eclipseState->getEclipseGrid()->getCartesianSize(), + 0)); + return eclWriter; +} + + +BOOST_AUTO_TEST_CASE(EclipseWriteRestartWellInfo) +{ + std::string eclipse_data_filename = "testBlackoilState3.DATA"; + std::string eclipse_restart_filename = "TESTBLACKOILSTATE3.X0004"; + + + test_work_area_type * test_area = test_work_area_alloc("TEST_EclipseWriteNumWells"); + test_work_area_copy_file(test_area, eclipse_data_filename.c_str()); + + Opm::ParseMode parseMode; + Opm::DeckConstPtr deck = createDeck(eclipse_data_filename); + Opm::EclipseStatePtr eclipseState(new Opm::EclipseState(deck , parseMode)); + Opm::EclipseWriterPtr eclipseWriter = createEclipseWriter(deck, eclipseState, eclipse_data_filename); + + std::shared_ptr simTimer( new Opm::SimulatorTimer() ); + simTimer->init(eclipseState->getSchedule()->getTimeMap()); + + eclipseWriter->writeInit(*simTimer); + + std::shared_ptr wellState(new Opm::WellState()); + std::shared_ptr blackoilState = createBlackOilState(eclipseState->getEclipseGrid()); + wellState->init(0, *blackoilState); + + int countTimeStep = eclipseState->getSchedule()->getTimeMap()->numTimesteps(); + + for(int timestep=0; timestep <= countTimeStep; ++timestep){ + simTimer->setCurrentStepNum(timestep); + eclipseWriter->writeTimeStep(*simTimer, *blackoilState, *wellState, false); + } + + verifyWellState(eclipse_restart_filename, eclipseState->getEclipseGrid(), eclipseState->getSchedule()); + + test_work_area_free(test_area); +}