OilfieldToolsNet/opm-core
4
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Deleted all files in opm/core/io

Browse Source 
This removes OutputWriter and eclipse, vtk, vag as this
functionality has moved to opm-output.
This commit is contained in:
Kjell W. Kongsvik 2016-03-17 14:15:50 +01:00
parent 03d5a5ba25
commit 8e68a4d815
20 changed files with 0 additions and 4585 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

24
CMakeLists_files.cmake
View File

@ -45,14 +45,6 @@ list (APPEND MAIN_SOURCE_FILES
opm/core/grid/cpgpreprocess/uniquepoints.c
opm/core/grid/grid.c
opm/core/grid/grid_equal.cpp
opm/core/io/OutputWriter.cpp
opm/core/io/eclipse/EclipseGridInspector.cpp
opm/core/io/eclipse/EclipseReader.cpp
opm/core/io/eclipse/EclipseWriteRFTHandler.cpp
opm/core/io/eclipse/EclipseWriter.cpp
opm/core/io/eclipse/writeECLData.cpp
opm/core/io/vag/vag.cpp
opm/core/io/vtk/writeVtkData.cpp
opm/core/linalg/LinearSolverFactory.cpp
opm/core/linalg/LinearSolverInterface.cpp
opm/core/linalg/LinearSolverIstl.cpp
@ -151,10 +143,6 @@ list (APPEND MAIN_SOURCE_FILES
# originally generated with the command:
# find tests -name '*.cpp' -a ! -wholename '*/not-unit/*' -printf '\t%p\n' | sort
list (APPEND TEST_SOURCE_FILES
tests/test_writenumwells.cpp
tests/test_writeReadRestartFile.cpp
tests/test_EclipseWriter.cpp
tests/test_EclipseWriteRFTHandler.cpp
tests/test_compressedpropertyaccess.cpp
tests/test_dgbasis.cpp
tests/test_cartgrid.cpp
@ -220,7 +208,6 @@ list (APPEND TEST_DATA_FILES
tests/satfuncEPS_D.DATA
tests/testBlackoilState1.DATA
tests/testBlackoilState2.DATA
tests/testBlackoilState3.DATA
tests/testPinch1.DATA
tests/wells_manager_data.data
tests/wells_manager_data_expanded.data
@ -294,17 +281,6 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/grid/cpgpreprocess/geometry.h
opm/core/grid/cpgpreprocess/preprocess.h
opm/core/grid/cpgpreprocess/uniquepoints.h
opm/core/io/OutputWriter.hpp
opm/core/io/eclipse/CornerpointChopper.hpp
opm/core/io/eclipse/EclipseGridInspector.hpp
opm/core/io/eclipse/EclipseIOUtil.hpp
opm/core/io/eclipse/EclipseReader.hpp
opm/core/io/eclipse/EclipseUnits.hpp
opm/core/io/eclipse/EclipseWriteRFTHandler.hpp
opm/core/io/eclipse/EclipseWriter.hpp
opm/core/io/eclipse/writeECLData.hpp
opm/core/io/vag/vag.hpp
opm/core/io/vtk/writeVtkData.hpp
opm/core/linalg/LinearSolverFactory.hpp
opm/core/linalg/LinearSolverInterface.hpp
opm/core/linalg/LinearSolverIstl.hpp

103
opm/core/io/OutputWriter.cpp
View File

@ -1,103 +0,0 @@
#include "OutputWriter.hpp"
#include <opm/core/grid.h>
#include <opm/core/io/eclipse/EclipseWriter.hpp>
#include <opm/core/utility/parameters/Parameter.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <forward_list>
#include <map>
#include <memory> // 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 <unique_ptr <OutputWriter> > writers_t;
typedef writers_t::iterator it_t;
typedef unique_ptr <writers_t> 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 SimulationDataContainer& 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 <typename Format> unique_ptr <OutputWriter>
create (const ParameterGroup& params,
std::shared_ptr <const EclipseState> eclipseState,
const Opm::PhaseUsage &phaseUsage,
std::shared_ptr <const UnstructuredGrid> grid) {
return unique_ptr <OutputWriter> (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 char*, unique_ptr <OutputWriter> (*)(
const ParameterGroup&,
std::shared_ptr <const EclipseState> eclipseState,
const Opm::PhaseUsage &phaseUsage,
std::shared_ptr <const UnstructuredGrid>)> map_t;
map_t FORMATS = {
{ "output_ecl", &create <EclipseWriter> },
};
} // anonymous namespace
unique_ptr <OutputWriter>
OutputWriter::create (const ParameterGroup& params,
std::shared_ptr <const EclipseState> eclipseState,
const Opm::PhaseUsage &phaseUsage,
std::shared_ptr <const UnstructuredGrid> 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 <bool> (name, false)) {
list->push_front (it->second (params, eclipseState, phaseUsage, grid));
}
}
// create a multiplexer from the list of formats we found
return unique_ptr <OutputWriter> (new MultiWriter (std::move (list)));
}

118
opm/core/io/OutputWriter.hpp
View File

@ -1,118 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_OUTPUT_WRITER_HPP
#define OPM_OUTPUT_WRITER_HPP
#include <memory> // unique_ptr, shared_ptr
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
struct UnstructuredGrid;
namespace Opm {
// forward declaration
class EclipseState;
namespace parameter { class ParameterGroup; }
class SimulationDataContainer;
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 <const Deck> (
* params.get <string> ("deck_filename"));
*
* std::unique_ptr <OutputWriter> 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 SimulationDataContainer& 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 <OutputWriter>
create (const parameter::ParameterGroup& params,
std::shared_ptr <const EclipseState> eclipseState,
const Opm::PhaseUsage &phaseUsage,
std::shared_ptr <const UnstructuredGrid> grid);
};
} // namespace Opm
#endif /* OPM_OUTPUT_WRITER_HPP */

457
opm/core/io/eclipse/CornerpointChopper.hpp
View File

@ -1,457 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
#define OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Parser/ParseContext.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <opm/parser/eclipse/Deck/DeckItem.hpp>
#include <iostream>
#include <fstream>
#include <string>
#include <vector>
#include <stdexcept>
#include <memory>
namespace Opm
{
class CornerPointChopper
{
public:
CornerPointChopper(const std::string& file)
{
Opm::ParseContext parseContext;
Opm::ParserPtr parser(new Opm::Parser());
deck_ = parser->parseFile(file , parseContext);
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<double>& 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<double, double> zLimits() const
{
return std::make_pair(botmax_, topmin_);
}
const std::pair<double, double> 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<double>& 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<double>& 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<Opm::UnitSystem> metricUnits_;
double botmax_;
double topmin_;
double abszmin_;
double abszmax_;
std::vector<double> new_COORD_;
std::vector<double> new_ZCORN_;
std::vector<int> new_ACTNUM_;
std::vector<double> new_PORO_;
std::vector<double> new_NTG_;
std::vector<double> new_SWCR_;
std::vector<double> new_SOWCR_;
std::vector<double> new_PERMX_;
std::vector<double> new_PERMY_;
std::vector<double> new_PERMZ_;
std::vector<int> new_SATNUM_;
int dims_[3];
int new_dims_[3];
std::vector<int> new_to_old_cell_;
void addDoubleKeyword_(Opm::DeckPtr subDeck,
const std::string& keywordName,
const std::string& dimensionString,
const std::vector<double>& 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<const Dimension> 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<int>& 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 <typename T>
void outputField(std::ostream& os,
const std::vector<T>& field,
const std::string& keyword,
const typename std::vector<T>::size_type nl = 20)
{
if (field.empty()) return;
os << keyword << '\n';
typedef typename std::vector<T>::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 <typename T>
void filterField(const std::vector<T>& field,
std::vector<T>& 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<double>& output_field)
{
if (deck_->hasKeyword(keyword)) {
const std::vector<double>& field = deck_->getKeyword(keyword).getRawDoubleData();
filterField(field, output_field);
}
}
void filterIntegerField(const std::string& keyword, std::vector<int>& output_field)
{
if (deck_->hasKeyword(keyword)) {
const std::vector<int>& field = deck_->getKeyword(keyword).getIntData();
filterField(field, output_field);
}
}
};
}
#endif // OPM_CORNERPOINTCHOPPER_HEADER_INCLUDED

350
opm/core/io/eclipse/EclipseGridInspector.cpp
View File

@ -1,350 +0,0 @@
//===========================================================================
//
// File: EclipseGridInspector.C
//
// Created: Mon Jun 2 12:17:51 2008
//
// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <opm/core/io/eclipse/EclipseGridInspector.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Deck/DeckItem.hpp>
#include <opm/parser/eclipse/Deck/DeckKeyword.hpp>
#include <opm/parser/eclipse/Deck/DeckRecord.hpp>
#include <stdexcept>
#include <numeric>
#include <cmath>
#include <cfloat>
#include <algorithm>
#include <array>
#include <iostream>
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<double,double> with x-dip in first component and y-dip in second.
*/
std::pair<double,double> EclipseGridInspector::cellDips(int i, int j, int k) const
{
checkLogicalCoords(i, j, k);
const std::vector<double>& 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<double>& 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<double, 8> 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<double, 6> 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<double,double> EclipseGridInspector::cellDips(int cell_idx) const
{
std::array<int, 3> idxs = cellIdxToLogicalCoords(cell_idx);
return cellDips(idxs[0], idxs[1], idxs[2]);
}
std::array<int, 3> 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<int, 3> a = {{i-1, j-1, k-1}};
return a; //std::array<int, 3> {{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<double>& 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<double>& 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<int, 3> 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<double, 6> 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<double> coord = deck_->getKeyword("COORD").getSIDoubleData();
std::vector<double> 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<double, 6> gridlimits = {{ xmin, xmax, ymin, ymax,
*min_element(zcorn.begin(), zcorn.end()),
*max_element(zcorn.begin(), zcorn.end()) }};
return gridlimits;
}
std::array<int, 3> EclipseGridInspector::gridSize() const
{
std::array<int, 3> retval = {{ logical_gridsize_[0],
logical_gridsize_[1],
logical_gridsize_[2] }};
return retval;
}
std::array<double, 8> EclipseGridInspector::cellZvals(int i, int j, int k) const
{
// Get the zcorn field.
const std::vector<double>& 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<double, 8> 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

103
opm/core/io/eclipse/EclipseGridInspector.hpp
View File

@ -1,103 +0,0 @@
//===========================================================================
//
// File: EclipseGridInspector.h
//
// Created: Mon Jun 2 09:46:08 2008
//
// Author: Atgeirr F Rasmussen <atgeirr@sintef.no>
//
// $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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLIPSEGRIDINSPECTOR_HEADER
#define OPM_ECLIPSEGRIDINSPECTOR_HEADER
#include <vector>
#include <array>
#include <opm/parser/eclipse/Deck/Deck.hpp>
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 <atgeirr@sintef.no>
@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<double,double> cellDips(int i, int j, int k) const;
std::pair<double,double> cellDips(int cell_idx) const;
// Convert global cell index to logical ijk-coordinates
std::array<int, 3> 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<double, 6> getGridLimits() const;
/// Returns the extent of the logical cartesian grid
/// as number of cells in the (i, j, k) directions.
std::array<int, 3> 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<double, 8> 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

34
opm/core/io/eclipse/EclipseIOUtil.hpp
View File

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

247
opm/core/io/eclipse/EclipseReader.cpp
View File

@ -1,247 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "EclipseReader.hpp"
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/grid/GridHelpers.hpp>
#include <opm/core/io/eclipse/EclipseIOUtil.hpp>
#include <opm/parser/eclipse/EclipseState/IOConfig/IOConfig.hpp>
#include <opm/parser/eclipse/EclipseState/InitConfig/InitConfig.hpp>
#include <opm/parser/eclipse/EclipseState/SimulationConfig/SimulationConfig.hpp>
#include <opm/parser/eclipse/Units/Dimension.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <algorithm>
#include <ert/ecl/ecl_file.h>
namespace Opm
{
static void restoreTemperatureData(const ecl_file_type* file,
EclipseStateConstPtr eclipse_state,
int numcells,
SimulationDataContainer& 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,
SimulationDataContainer& 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");
}
}
static void restoreSaturation(const ecl_file_type* file_type,
const PhaseUsage& phaseUsage,
int numcells,
SimulationDataContainer& 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<double> 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<double> sgas_datavec(&sgas_data[0], &sgas_data[numcells]);
EclipseIOUtil::addToStripedData(sgas_datavec, simulator_state.saturation(), phaseUsage.phase_pos[BlackoilPhases::Vapour], phaseUsage.num_phases);
} else {
std::string error_str = "Restart file is missing SGAS data!\n";
throw std::runtime_error(error_str);
}
}
}
static void restoreRSandRV(const ecl_file_type* file_type,
SimulationConfigConstPtr sim_config,
int numcells,
SimulationDataContainer& simulator_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);
auto& rs = simulator_state.getCellData( BlackoilState::GASOILRATIO );
for (int i = 0; i < ecl_kw_get_size( rs_kw ); i++) {
rs[i] = rs_data[i];
}
} 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);
auto& rv = simulator_state.getCellData( BlackoilState::RV );
for (int i = 0; i < ecl_kw_get_size( rv_kw ); i++) {
rv[i] = rv_data[i];
}
} else {
throw std::runtime_error("Restart file is missing RV data!\n");
}
}
}
static void restoreSOLUTION(const std::string& restart_filename,
int reportstep,
bool unified,
EclipseStateConstPtr eclipseState,
int numcells,
const PhaseUsage& phaseUsage,
SimulationDataContainer& 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);
if (simulator_state.hasCellData( BlackoilState::RV )) {
SimulationConfigConstPtr sim_config = eclipseState->getSimulationConfig();
restoreRSandRV(file_type, sim_config, numcells, simulator_state );
}
} 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);
}
}
static 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,
SimulationDataContainer& 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

39
opm/core/io/eclipse/EclipseReader.hpp
View File

@ -1,39 +0,0 @@
#ifndef ECLIPSEREADER_HPP
#define ECLIPSEREADER_HPP
#include <string>
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
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<double> objects
///
class SimulationDataContainer;
void init_from_restart_file(EclipseStateConstPtr eclipse_state,
int numcells,
const PhaseUsage& pu,
SimulationDataContainer& simulator_state,
WellState& wellstate);
}
#endif // ECLIPSEREADER_HPP

63
opm/core/io/eclipse/EclipseUnits.hpp
View File

@ -1,63 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#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

144
opm/core/io/eclipse/EclipseWriteRFTHandler.cpp
View File

@ -1,144 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include <vector>
#include <opm/core/io/eclipse/EclipseWriteRFTHandler.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/core/utility/miscUtilities.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/CompletionSet.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/WellSet.hpp>
#include <opm/parser/eclipse/EclipseState/Grid/EclipseGrid.hpp>
#include <ert/ecl/ecl_rft_node.h>
#include <ert/ecl/ecl_rft_file.h>
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<WellConstPtr>& wells,
EclipseGridConstPtr eclipseGrid,
std::vector<double>& pressure,
std::vector<double>& swat,
std::vector<double>& sgas) {
std::vector<ecl_rft_node_type *> rft_nodes;
for (std::vector<WellConstPtr>::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<double>& pressure,
const std::vector<double>& swat,
const std::vector<double>& 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

77
opm/core/io/eclipse/EclipseWriteRFTHandler.hpp
View File

@ -1,77 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLIPSE_WRITE_RFT_HANDLER_HPP
#define OPM_ECLIPSE_WRITE_RFT_HANDLER_HPP
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <ert/ecl/ecl_rft_node.h>
#include <ert/ecl/ecl_util.h>
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<std::shared_ptr< const Well >>& wells,
std::shared_ptr< const EclipseGrid > eclipseGrid,
std::vector<double>& pressure,
std::vector<double>& swat,
std::vector<double>& 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<double>& pressure,
const std::vector<double>& swat,
const std::vector<double>& sgas);
void initGlobalToActiveIndex(const int * compressedToCartesianCellIdx, size_t numCells, size_t cartesianSize);
std::vector<int> globalToActiveIndex_;
};
}//namespace EclipseWriterDetails
}//namespace Opm
#endif //OPM_ECLIPSE_WRITE_RFT_HANDLER_HPP

1508
opm/core/io/eclipse/EclipseWriter.cpp
View File

File diff suppressed because it is too large Load Diff

139
opm/core/io/eclipse/EclipseWriter.hpp
View File

@ -1,139 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_ECLIPSE_WRITER_HPP
#define OPM_ECLIPSE_WRITER_HPP
#include <opm/core/io/OutputWriter.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/wells.h> // WellType
#include <opm/core/simulator/SimulatorTimerInterface.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/Units/UnitSystem.hpp>
#include <ert/ecl/ecl_util.h>
#include <string>
#include <vector>
#include <array>
#include <memory>
struct UnstructuredGrid;
namespace Opm {
// forward declarations
namespace EclipseWriterDetails {
class Summary;
}
class SimulationDataContainer;
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 SimulationDataContainer& 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<int, 3> 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<EclipseWriterDetails::Summary> summary_;
void init(const parameter::ParameterGroup& params);
};
typedef std::shared_ptr<EclipseWriter> EclipseWriterPtr;
typedef std::shared_ptr<const EclipseWriter> EclipseWriterConstPtr;
} // namespace Opm
#endif // OPM_ECLIPSE_WRITER_HPP

193
opm/core/io/eclipse/writeECLData.cpp
View File

@ -1,193 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif
#include <opm/core/grid.h>
#include <opm/core/io/eclipse/writeECLData.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/common/ErrorMacros.hpp>
#include <vector>
#ifdef HAVE_ERT // This one goes almost to the bottom of the file
#include <ert/ecl/ecl_grid.h>
#include <ert/ecl/ecl_util.h>
#include <ert/ecl/ecl_rst_file.h>
namespace Opm
{
static ecl_kw_type * ecl_kw_wrapper( const UnstructuredGrid& grid,
const std::string& kw_name ,
const std::vector<double> * 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<double>::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

46
opm/core/io/eclipse/writeECLData.hpp
View File

@ -1,46 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WRITEECLDATA_HEADER_INCLUDED
#define OPM_WRITEECLDATA_HEADER_INCLUDED
#include <opm/core/utility/DataMap.hpp>
#include <boost/date_time/posix_time/posix_time_types.hpp>
#include <string>
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

408
opm/core/io/vag/vag.cpp
View File

@ -1,408 +0,0 @@
/*===========================================================================
//
// File: vag.cpp
//
// Created: 2012-06-08 15:45:53+0200
//
// Authors: Knut-Andreas Lie <Knut-Andreas.Lie@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Xavier Raynaud <Xavier.Raynaud@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/io/vag/vag.hpp>
#include <opm/core/grid/cornerpoint_grid.h>
#include <iostream>
#include <fstream>
#include <string>
#include <stdlib.h>
#include <cmath>
#include <cassert>
#include <set>
#include <vector>
#include <map>
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 <<number << endl;
pos_struct.pos[i+1]=pos_struct.pos[i]+number;
for(int j=0;j< number;++j){
int value;
is >> 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<int> > 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<int> 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<int>::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<int> > edges;
std::vector< std::vector< std::set<int> > > 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<int> spair;
spair.insert(node1);
spair.insert(node2);
edges.insert(spair);
faces_spares[i].push_back(spair);
}
// add end segment
{
std::set<int> 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<std::set<int>, int> edge_map;
std::set< std::set<int> >::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> , int >(*it,k));
k=k+1;
std::set<int>::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;j<int(faces_spares[i].size());++j){
int edge_num=edge_map[faces_spares[i][j]];
vag_grid.faces_to_edges.value.push_back(edge_num+1);
}
}
// vag_grid.edges(0);//not known
// vag_grid.faces_to_edges// not known
// material // can not be extracted from the grid
}
void writeVagFormat(std::ostream& os,Opm::VAG& vag_grid){
using namespace std;
os << "File in the Vag grid format\n";
os << "Number of vertices " ;
os << vag_grid.number_of_vertices << endl;;
os <<"Number of control volume ";
os << vag_grid.number_of_volumes << endl;
os <<"Number of faces " ;
os << vag_grid.number_of_faces << endl;
os <<"Number of edges " ;
os << vag_grid.number_of_edges << endl;
os <<"Vertices " << vag_grid.vertices.size() << endl;
writeVector(os, vag_grid.vertices,3);
os << "Volumes->faces " << 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);
*/
}
}

165
opm/core/io/vag/vag.hpp
View File

@ -1,165 +0,0 @@
/*===========================================================================
//
// File: vag.hpp
//
// Created: 2012-06-08 15:46:23+0200
//
// Authors: Knut-Andreas Lie <Knut-Andreas.Lie@sintef.no>
// Halvor M. Nilsen <HalvorMoll.Nilsen@sintef.no>
// Atgeirr F. Rasmussen <atgeirr@sintef.no>
// Xavier Raynaud <Xavier.Raynaud@sintef.no>
// Bård Skaflestad <Bard.Skaflestad@sintef.no>
//
//==========================================================================*/
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_VAG_HPP_HEADER
#define OPM_VAG_HPP_HEADER
#include <opm/core/grid.h>
#include <istream>
#include <ostream>
#include <string>
#include <vector>
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<int> pos;
std::vector<int> 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<double> 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<int> edges;
/** The two neigbours of the face i is faces_to_volumes[2*i:2*i+1] */
std::vector<int> 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<double> 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 <typename T>
void readVector(std::istream& is,std::vector<T>& 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 <typename T>
void writeVector(std::ostream& os,std::vector<T>& vec,int n){
typedef typename std::vector<T>::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 */

319
opm/core/io/vtk/writeVtkData.cpp
View File

@ -1,319 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include "config.h"
#include <opm/core/io/vtk/writeVtkData.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/grid.h>
#include <opm/common/ErrorMacros.hpp>
#include <boost/lexical_cast.hpp>
#include <set>
#include <cmath>
#include <algorithm>
#include <iostream>
#include <iterator>
#include <vector>
namespace Opm
{
void writeVtkData(const std::array<int, 3>& dims,
const std::array<double, 3>& 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<double>& 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<std::string, std::string> 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_ << "</" << name_ << ">\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 << "<?xml version=\"1.0\"?>\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<std::string>(num_pts);
pm["NumberOfCells"] = boost::lexical_cast<std::string>(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<int> 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<int> 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<int>(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<int> 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<int>(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<double>& field = *(dit->second);
const int num_comps = field.size()/grid.number_of_cells;
pm["NumberOfComponents"] = boost::lexical_cast<std::string>(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<double>::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

48
opm/core/io/vtk/writeVtkData.hpp
View File

@ -1,48 +0,0 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_WRITEVTKDATA_HEADER_INCLUDED
#define OPM_WRITEVTKDATA_HEADER_INCLUDED
#include <string>
#include <map>
#include <vector>
#include <array>
#include <iosfwd>
#include <opm/core/utility/DataMap.hpp>
struct UnstructuredGrid;
namespace Opm
{
/// Vtk output for cartesian grids.
void writeVtkData(const std::array<int, 3>& dims,
const std::array<double, 3>& 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