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

Merge pull request #67 from bska/master

Browse Source 
Merge ERT Branch
This commit is contained in:
Atgeirr Flø Rasmussen 2012-10-12 10:38:51 -07:00
commit 81f8a9895d
17 changed files with 1176 additions and 289 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

13
Makefile.am
View File

@ -12,6 +12,7 @@ lib_LTLIBRARIES = lib/libopmcore.la
# Build-time flags needed to build libopmcore.la # Build-time flags needed to build libopmcore.la
AM_CPPFLAGS = \ AM_CPPFLAGS = \
$(ERT_CPPFLAGS) \
$(OPM_BOOST_CPPFLAGS) \ $(OPM_BOOST_CPPFLAGS) \
${SUPERLU_CPPFLAGS} ${SUPERLU_CPPFLAGS}
@ -22,6 +23,7 @@ ${SUPERLU_CPPFLAGS}
lib_libopmcore_la_LDFLAGS = \ lib_libopmcore_la_LDFLAGS = \
-R $(OPM_BOOST_LIBDIR) \ -R $(OPM_BOOST_LIBDIR) \
$(OPM_BOOST_LDFLAGS) \ $(OPM_BOOST_LDFLAGS) \
$(ERT_LDFLAGS)
${SUPERLU_LDFLAGS} ${SUPERLU_LDFLAGS}
lib_libopmcore_la_LIBADD = \ lib_libopmcore_la_LIBADD = \
@ -29,6 +31,7 @@ $(BOOST_FILESYSTEM_LIB) \
$(BOOST_SYSTEM_LIB) \ $(BOOST_SYSTEM_LIB) \
$(BOOST_DATE_TIME_LIB) \ $(BOOST_DATE_TIME_LIB) \
$(BOOST_UNIT_TEST_FRAMEWORK_LIB) \ $(BOOST_UNIT_TEST_FRAMEWORK_LIB) \
$(ERT_LIBS) \
$(LAPACK_LIBS) ${SUPERLU_LIBS} $(BLAS_LIBS) $(LIBS) $(LAPACK_LIBS) ${SUPERLU_LIBS} $(BLAS_LIBS) $(LIBS)
# ---------------------------------------------------------------------- # ----------------------------------------------------------------------
@ -242,6 +245,7 @@ opm/core/transport/spu_implicit.h \
opm/core/transport/transport_source.h \ opm/core/transport/transport_source.h \
opm/core/utility/Average.hpp \ opm/core/utility/Average.hpp \
opm/core/utility/ColumnExtract.hpp \ opm/core/utility/ColumnExtract.hpp \
opm/core/utility/DataMap.hpp \
opm/core/utility/ErrorMacros.hpp \ opm/core/utility/ErrorMacros.hpp \
opm/core/utility/Factory.hpp \ opm/core/utility/Factory.hpp \
opm/core/utility/MonotCubicInterpolator.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \
@ -292,6 +296,15 @@ opm/core/linalg/LinearSolverUmfpack.hpp
endif endif
if HAVE_ERT
lib_libopmcore_la_SOURCES += \
opm/core/utility/writeECLData.cpp
nobase_include_HEADERS += \
opm/core/utility/writeECLData.hpp
endif
if DUNE_ISTL if DUNE_ISTL
lib_libopmcore_la_SOURCES += \ lib_libopmcore_la_SOURCES += \
opm/core/linalg/LinearSolverIstl.cpp opm/core/linalg/LinearSolverIstl.cpp

2
configure.ac
View File

@ -39,6 +39,8 @@ OPM_CORE_CHECKS
OPM_DYNLINK_BOOST_TEST OPM_DYNLINK_BOOST_TEST
ERT
dnl Substitute Autoconf's abs_*dir variables into the Makefiles for the dnl Substitute Autoconf's abs_*dir variables into the Makefiles for the
dnl benefit of external code that uses these variables to derive dnl benefit of external code that uses these variables to derive
dnl locations (e.g., Dune's DUNE_CHECK_MODULES macro). Automakes prior dnl locations (e.g., Dune's DUNE_CHECK_MODULES macro). Automakes prior

3
examples/Makefile.am
View File

@ -1,6 +1,9 @@
# Build-time flags needed to form example programs # Build-time flags needed to form example programs
ERT_INCLUDE_PATH = $(ERT_ROOT)/include
AM_CPPFLAGS = \ AM_CPPFLAGS = \
-I$(top_srcdir) \ -I$(top_srcdir) \
-I$(ERT_INCLUDE_PATH) \
$(OPM_BOOST_CPPFLAGS) $(OPM_BOOST_CPPFLAGS)
# All targets link to the library # All targets link to the library

14
examples/spu_2p.cpp
View File

@ -94,14 +94,19 @@
#include <numeric> #include <numeric>
#ifdef HAVE_ERT
#include <opm/core/utility/writeECLData.hpp>
#endif
static void outputState(const UnstructuredGrid& grid, static void outputState(const UnstructuredGrid& grid,
const Opm::TwophaseState& state, const Opm::TwophaseState& state,
const int step, const Opm::SimulatorTimer& simtimer,
const std::string& output_dir) const std::string& output_dir)
{ {
// Write data in VTK format. // Write data in VTK format.
int step = simtimer.currentStepNum();
std::ostringstream vtkfilename; std::ostringstream vtkfilename;
vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu";
std::ofstream vtkfile(vtkfilename.str().c_str()); std::ofstream vtkfile(vtkfilename.str().c_str());
@ -115,6 +120,9 @@ static void outputState(const UnstructuredGrid& grid,
Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity);
dm["velocity"] = &cell_velocity; dm["velocity"] = &cell_velocity;
Opm::writeVtkData(grid, dm, vtkfile); Opm::writeVtkData(grid, dm, vtkfile);
#ifdef HAVE_ERT
Opm::writeECLData(grid , dm , simtimer , output_dir , "OPM" );
#endif
// Write data (not grid) in Matlab format // Write data (not grid) in Matlab format
for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) { for (Opm::DataMap::const_iterator it = dm.begin(); it != dm.end(); ++it) {
@ -530,7 +538,7 @@ main(int argc, char** argv)
// Report timestep and (optionally) write state to disk. // Report timestep and (optionally) write state to disk.
simtimer.report(std::cout); simtimer.report(std::cout);
if (output && (simtimer.currentStepNum() % output_interval == 0)) { if (output && (simtimer.currentStepNum() % output_interval == 0)) {
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); outputState(*grid->c_grid(), state, simtimer , output_dir);
} }
// Solve pressure. // Solve pressure.
@ -699,7 +707,7 @@ main(int argc, char** argv)
<< "\n Transport time: " << ttime << std::endl; << "\n Transport time: " << ttime << std::endl;
if (output) { if (output) {
outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); outputState(*grid->c_grid(), state, simtimer, output_dir);
outputWaterCut(watercut, output_dir); outputWaterCut(watercut, output_dir);
if (wells->c_wells()) { if (wells->c_wells()) {
outputWellReport(wellreport, output_dir); outputWellReport(wellreport, output_dir);

67
m4/ert.m4 Normal file
View File

@ -0,0 +1,67 @@
AC_DEFUN([_ERT_SOURCE_TEXT],
[
AC_LANG_PROGRAM(
[[
#include <stddef.h>
#include <ecl_util.h>
]],dnl
[[
int sz;
sz = ecl_util_get_sizeof_ctype(ECL_INT_TYPE);
]])[]dnl
])[]dnl
# ----------------------------------------------------------------------
AC_DEFUN([ERT],
[
AC_ARG_WITH([ert],
[AS_HELP_STRING([--with-ert=<root>], [Use ERT libraries])],
[], [with_ert=no])
use_ert=no
AS_IF([test x"${with_ert}" != x"no"],
[
_ert_LDFLAGS_SAVE="${LDFLAGS}"
_ert_LIBS_SAVE="${LIBS}"
_ert_CPPFLAGS_SAVE="${CPPFLAGS}"
_ert_CFLAGS_SAVE="${CFLAGS}"
ERT_CPPFLAGS=
ERT_LDFLAGS=
ERT_LIBS="-lecl -lgeometry -lert_util -lpthread -lz -lgomp"
AS_IF([test x"${with_ert}" != x"yes"],
[ERT_LDFLAGS="-L${with_ert}/lib"
ERT_CPPFLAGS="-I${with_ert}/include"], [:])[]dnl
CFLAGS="-std=gnu99"
CPPFLAGS="${ERT_CPPFLAGS} ${CPPFLAGS}"
LDFLAGS="${ERT_LDFLAGS} ${LDFLAGS}"
LIBS="${ERT_LIBS} ${LIBS}"
AC_LINK_IFELSE([_ERT_SOURCE_TEXT],
[use_ert=yes], [use_ert=no])
LIBS="${_ert_LIBS_SAVE}"
CPPFLAGS="${_ert_CPPFLAGS_SAVE}"
LDFLAGS="${_ert_LDFLAGS_SAVE}"
CFLAGS="${_ert_CFLAGS_SAVE}"
AS_IF([test x"${use_ert}" == x"yes"],
[AC_SUBST([ERT_CPPFLAGS])
AC_SUBST([ERT_LDFLAGS])
AC_SUBST([ERT_LIBS])
AC_DEFINE([HAVE_ERT], [1],
[Are the ERT libraries available for reading and writing ECLIPSE files.])],dnl
[:])[]dnl
], [:])[]dnl
AM_CONDITIONAL([HAVE_ERT], [test x"${use_ert}" != x"no"])
# AC_MSG_ERROR(
# [**** ERT_CPPFLAGS = ${ERT_CPPFLAGS} ****
# **** ERT_LDFLAGS = ${ERT_LDFLAGS} ****
# **** ERT_LIBS = ${ERT_LIBS} ****
# ])
])

26
opm/core/GridManager.cpp
View File

@ -37,7 +37,8 @@ namespace Opm
{ {
// We accept two different ways to specify the grid. // We accept two different ways to specify the grid.
// 1. Corner point format. // 1. Corner point format.
// Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally ACTNUM. // Requires ZCORN, COORDS, DIMENS or SPECGRID, optionally
// ACTNUM, optionally MAPAXES.
// For this format, we will verify that DXV, DYV, DZV, // For this format, we will verify that DXV, DYV, DZV,
// DEPTHZ and TOPS are not present. // DEPTHZ and TOPS are not present.
// 2. Tensor grid format. // 2. Tensor grid format.
@ -119,29 +120,8 @@ namespace Opm
void GridManager::initFromDeckCornerpoint(const Opm::EclipseGridParser& deck) void GridManager::initFromDeckCornerpoint(const Opm::EclipseGridParser& deck)
{ {
// Extract data from deck. // Extract data from deck.
const std::vector<double>& zcorn = deck.getFloatingPointValue("ZCORN");
const std::vector<double>& coord = deck.getFloatingPointValue("COORD");
const int* actnum = 0;
if (deck.hasField("ACTNUM")) {
actnum = &(deck.getIntegerValue("ACTNUM")[0]);
}
std::vector<int> dims;
if (deck.hasField("DIMENS")) {
dims = deck.getIntegerValue("DIMENS");
} else if (deck.hasField("SPECGRID")) {
dims = deck.getSPECGRID().dimensions;
} else {
THROW("Deck must have either DIMENS or SPECGRID.");
}
// Collect in input struct for preprocessing. // Collect in input struct for preprocessing.
struct grdecl grdecl; struct grdecl grdecl = deck.get_grdecl();
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = actnum;
grdecl.dims[0] = dims[0];
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];
// Process grid. // Process grid.
ug_ = create_grid_cornerpoint(&grdecl, 0.0); ug_ = create_grid_cornerpoint(&grdecl, 0.0);

339
opm/core/eclipse/EclipseGridParser.cpp
View File

@ -47,8 +47,19 @@
#include <opm/core/eclipse/EclipseGridParserHelpers.hpp> #include <opm/core/eclipse/EclipseGridParserHelpers.hpp>
#include <opm/core/eclipse/SpecialEclipseFields.hpp> #include <opm/core/eclipse/SpecialEclipseFields.hpp>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <boost/filesystem.hpp>
#include <opm/core/utility/Units.hpp> #include <opm/core/utility/Units.hpp>
#include <opm/core/grid/cpgpreprocess/preprocess.h>
#include <boost/filesystem.hpp>
#ifdef HAVE_ERT
#include <fortio.h>
#include <ecl_grid.h>
#include <ecl_kw_magic.h>
#include <ecl_kw.h>
#include <ecl_util.h>
#include <ecl_init_file.h>
#include <ecl_file.h>
#endif
using namespace std; using namespace std;
@ -86,7 +97,7 @@ namespace EclipseKeywords
string("MULTPV"), string("PRESSURE"), string("SGAS"), string("MULTPV"), string("PRESSURE"), string("SGAS"),
string("SWAT"), string("SOIL"), string("RS"), string("SWAT"), string("SOIL"), string("RS"),
string("DXV"), string("DYV"), string("DZV"), string("DXV"), string("DYV"), string("DZV"),
string("DEPTHZ"), string("TOPS") string("DEPTHZ"), string("TOPS"), string("MAPAXES")
}; };
const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]); const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
@ -114,7 +125,7 @@ namespace EclipseKeywords
const int num_special_fields = sizeof(special_fields) / sizeof(special_fields[0]); const int num_special_fields = sizeof(special_fields) / sizeof(special_fields[0]);
string ignore_with_data[] = string ignore_with_data[] =
{ string("MAPUNITS"), string("MAPAXES"), string("GRIDUNIT"), { string("MAPUNITS"), string("GRIDUNIT"),
string("NTG"), string("REGDIMS"), string("WELLDIMS"), string("NTG"), string("REGDIMS"), string("WELLDIMS"),
string("NSTACK"), string("SATNUM"), string("NSTACK"), string("SATNUM"),
string("RPTRST"), string("ROIP"), string("RWIP"), string("RPTRST"), string("ROIP"), string("RWIP"),
@ -141,6 +152,9 @@ namespace EclipseKeywords
string include_keywords[] = { string("INCLUDE") }; string include_keywords[] = { string("INCLUDE") };
const int num_include_keywords = sizeof(include_keywords) / sizeof(include_keywords[0]); const int num_include_keywords = sizeof(include_keywords) / sizeof(include_keywords[0]);
string import_keywords[] = { string("IMPORT") };
const int num_import_keywords = sizeof(import_keywords) / sizeof(import_keywords[0]);
} // namespace EclipseKeywords } // namespace EclipseKeywords
@ -154,6 +168,7 @@ namespace {
IgnoreWithData, IgnoreWithData,
IgnoreNoData, IgnoreNoData,
Include, Include,
Import,
Unknown Unknown
}; };
@ -174,6 +189,8 @@ namespace {
return IgnoreNoData; return IgnoreNoData;
} else if (count(include_keywords, include_keywords + num_include_keywords, keyword)) { } else if (count(include_keywords, include_keywords + num_include_keywords, keyword)) {
return Include; return Include;
} else if (count(import_keywords, import_keywords + num_import_keywords, keyword)) {
return Import;
} else { } else {
return Unknown; return Unknown;
} }
@ -488,6 +505,14 @@ void EclipseGridParser::readImpl(istream& is)
// is >> ignoreSlashLine; // is >> ignoreSlashLine;
break; break;
} }
case Import: {
string import_filename = readString(is);
if (!directory_.empty()) {
import_filename = directory_ + '/' + import_filename;
}
getNumericErtFields(import_filename);
break;
}
case Unknown: case Unknown:
default: default:
ignored_fields_.insert(keyword); ignored_fields_.insert(keyword);
@ -542,6 +567,9 @@ void EclipseGridParser::convertToSI()
do_convert = false; // Dimensionless keywords... do_convert = false; // Dimensionless keywords...
} else if (key == "PRESSURE") { } else if (key == "PRESSURE") {
unit = units_.pressure; unit = units_.pressure;
} else if (key == "MAPAXES") {
MESSAGE("Not applying units to MAPAXES yet!");
unit = 1.0;
} else { } else {
THROW("Units for field " << key << " not specified. Cannon convert to SI."); THROW("Units for field " << key << " not specified. Cannon convert to SI.");
} }
@ -800,5 +828,310 @@ void EclipseGridParser::computeUnits()
THROW("Unknown unit family " << unit_family); THROW("Unknown unit family " << unit_family);
} }
} }
struct grdecl EclipseGridParser::get_grdecl() const {
struct grdecl grdecl;
// Extract data from deck.
const std::vector<double>& zcorn = getFloatingPointValue("ZCORN");
const std::vector<double>& coord = getFloatingPointValue("COORD");
const int* actnum = NULL;
if (hasField("ACTNUM")) {
actnum = &(getIntegerValue("ACTNUM")[0]);
}
std::vector<int> dims;
if (hasField("DIMENS")) {
dims = getIntegerValue("DIMENS");
} else if (hasField("SPECGRID")) {
dims = getSPECGRID().dimensions;
} else {
THROW("Deck must have either DIMENS or SPECGRID.");
}
// Collect in input struct for preprocessing.
grdecl.zcorn = &zcorn[0];
grdecl.coord = &coord[0];
grdecl.actnum = actnum;
grdecl.dims[0] = dims[0];
grdecl.dims[1] = dims[1];
grdecl.dims[2] = dims[2];
if (hasField("MAPAXES")) {
const std::vector<double> &mapaxes = getFloatingPointValue("MAPAXES");
grdecl.mapaxes = &mapaxes[0];
} else
grdecl.mapaxes = NULL;
return grdecl;
}
#ifdef HAVE_ERT
/*
This function will create a ecl_kw instance filled with the data
from input argument @keyword. The ecl_kw will get it's own copy of
the data.
If the input type ecl_type == ECL_INT_TYPE the function will use the
getIntegerValue() method to get the keyword data, if ecl_type ==
ECL_DOUBLE_TYPE || ECL_FLOAT_TYPE the getFloatingPointValue()
function is invoked. If ecl_type == ECL_FLOAT_TYPE the data will be
converted to when inserting into the ecl_kw.
When the ecl_kw instance is no longer needed it should be discarded
with a call to ecl_kw_free( ).
If you are asking for a non-existent field the function will return NULL
*/
ecl_kw_type * EclipseGridParser::newEclKW(const std::string &keyword , ecl_type_enum ecl_type) const {
ecl_kw_type * ecl_kw = NULL;
if (hasField(keyword)) {
if (ecl_type == ECL_INT_TYPE) {
std::vector<int> data = getIntegerValue( keyword );
ecl_kw = ecl_kw_alloc( keyword.c_str() , data.size() , ecl_type );
ecl_kw_set_memcpy_data( ecl_kw , &data[0]);
} else {
std::vector<double> data = getFloatingPointValue( keyword );
if (ecl_type == ECL_DOUBLE_TYPE) {
ecl_kw = ecl_kw_alloc( keyword.c_str() , data.size() , ecl_type );
ecl_kw_set_memcpy_data( ecl_kw , &data[0]);
} else if (ecl_type == ECL_FLOAT_TYPE) {
ecl_kw = ecl_kw_alloc( keyword.c_str() , data.size() , ecl_type );
for (std::vector<double>::size_type i=0; i < data.size(); i++)
ecl_kw_iset_float( ecl_kw , i , data[i] );
}
}
}
return ecl_kw;
}
/**
This function will extract the COORD, ZCORN, ACTNUM and optionaly
MAPAXES keywords from the eclipse deck and create an ecl_grid
instance.
When you are finished working with the ecl_grid instance it should
be disposed with ecl_grid_free( ).
*/
ecl_grid_type * EclipseGridParser::newGrid( ) {
struct grdecl grdecl = get_grdecl();
ecl_kw_type * coord_kw = newEclKW( COORD_KW , ECL_FLOAT_TYPE );
ecl_kw_type * zcorn_kw = newEclKW( ZCORN_KW , ECL_FLOAT_TYPE );
ecl_kw_type * actnum_kw = newEclKW( ACTNUM_KW , ECL_INT_TYPE );
ecl_kw_type * mapaxes_kw = NULL;
ecl_grid_type * grid ;
if (grdecl.mapaxes != NULL)
mapaxes_kw = newEclKW( MAPAXES_KW , ECL_FLOAT_TYPE );
grid = ecl_grid_alloc_GRDECL_kw( grdecl.dims[0] , grdecl.dims[1] , grdecl.dims[2] , zcorn_kw , coord_kw , actnum_kw , mapaxes_kw );
ecl_kw_free( coord_kw );
ecl_kw_free( zcorn_kw );
ecl_kw_free( actnum_kw );
if (mapaxes_kw != NULL)
ecl_kw_free( mapaxes_kw );
return grid;
}
/**
This function will save an EGRID file based on the COORD, ZCORN,
ACTNUM and optionally MAPAXES keywords included in the deck.
This function creates the EGRID file without going through a
ecl_grid instance; this is obviously somewhat faster and less
memory demanding. Alternatively you can create a ecl_grid instance
and then subsequently store that grid as an EGRID file:
{
ecl_grid_type * grid = newGRID( );
ecl_grid_fwrite_EGRID( grid , filename );
ecl_grid_free( grid );
}
*/
void EclipseGridParser::saveEGRID( const std::string & filename) {
bool endian_flip = true;//ECL_ENDIAN_FLIP;
bool fmt_file = ecl_util_fmt_file( filename.c_str() );
struct grdecl grdecl = get_grdecl();
fortio_type * fortio = fortio_open_writer( filename.c_str() , fmt_file , endian_flip );
{
float * mapaxes = NULL;
if (grdecl.mapaxes != NULL) {
mapaxes = new float[6];
for (int i=0; i < 6; i++)
mapaxes[i]= grdecl.mapaxes[i];
}
ecl_grid_fwrite_EGRID_header( grdecl.dims , mapaxes , fortio );
if (grdecl.mapaxes != NULL)
delete[] mapaxes;
}
{
ecl_kw_type * coord_kw = newEclKW( COORD_KW , ECL_FLOAT_TYPE );
ecl_kw_type * zcorn_kw = newEclKW( ZCORN_KW , ECL_FLOAT_TYPE );
ecl_kw_type * actnum_kw = newEclKW( ACTNUM_KW , ECL_INT_TYPE );
ecl_kw_type * endgrid_kw = ecl_kw_alloc( ENDGRID_KW , 0 , ECL_INT_TYPE );
ecl_kw_fwrite( coord_kw , fortio );
ecl_kw_fwrite( zcorn_kw , fortio );
ecl_kw_fwrite( actnum_kw , fortio );
ecl_kw_fwrite( endgrid_kw , fortio );
ecl_kw_free( coord_kw );
ecl_kw_free( zcorn_kw );
ecl_kw_free( actnum_kw );
ecl_kw_free( endgrid_kw );
}
fortio_fclose( fortio );
}
/**
Will query the deck for keyword @kw; and save it to the @fortio
instance if the keyword can be found.
*/
void EclipseGridParser::save_kw( fortio_type * fortio , const std::string & kw , ecl_type_enum ecl_type) {
ecl_kw_type * ecl_kw = newEclKW( kw , ecl_type );
if (ecl_kw != NULL) {
ecl_kw_fwrite( ecl_kw , fortio );
ecl_kw_free( ecl_kw );
}
}
/**
Will save an ECLIPSE INIT file to @filename. Observe that the main
focus of this function is to store grid properties like PERMX and
PORO, various tabular properties like e.g. relperm tables and
thermodynamic properties are ignored.
*/
void EclipseGridParser::saveINIT( const std::string & filename , const ecl_grid_type * ecl_grid) {
int phases = ECL_OIL_PHASE + ECL_WATER_PHASE;
bool fmt_file = ecl_util_fmt_file( filename.c_str() );
bool endian_flip = true;//ECL_ENDIAN_FLIP;
fortio_type * fortio = fortio_open_writer( filename.c_str() , fmt_file , endian_flip );
{
ecl_kw_type * poro_kw = newEclKW( PORO_KW , ECL_FLOAT_TYPE );
time_t start_date;
{
tm td_tm = to_tm( start_date_ );
start_date = mktime( &td_tm );
}
ecl_init_file_fwrite_header( fortio , ecl_grid , poro_kw , phases , start_date );
ecl_kw_free( poro_kw );
}
/* This collection of keywords is somewhat arbitrary and random. */
save_kw( fortio , "PERMX" , ECL_FLOAT_TYPE);
save_kw( fortio , "PERMY" , ECL_FLOAT_TYPE);
save_kw( fortio , "PERMZ" , ECL_FLOAT_TYPE);
save_kw( fortio , "FIPNUM" , ECL_INT_TYPE);
save_kw( fortio , "SATNUM" , ECL_INT_TYPE);
save_kw( fortio , "EQLNUM" , ECL_INT_TYPE);
fortio_fclose( fortio );
}
/**
This is the main function used to save the state of the ECLIPSE
deck in ECLIPSE format. The function will save an INIT file and an
EGRID file.
The input arguments are the output directory to store files in, and
the basename to use for the files; the function will build up a
ECLIPSE standard filename internally.
*/
void EclipseGridParser::saveEGRID_INIT( const std::string& output_dir , const std::string& basename, bool fmt_file) {
ecl_grid_type * ecl_grid = newGrid();
char * egrid_file = ecl_util_alloc_filename( output_dir.c_str() , basename.c_str() , ECL_EGRID_FILE , fmt_file , 0);
char * init_file = ecl_util_alloc_filename( output_dir.c_str() , basename.c_str() , ECL_INIT_FILE , fmt_file , 0);
ecl_grid_fwrite_EGRID( ecl_grid , egrid_file );
saveINIT( init_file , ecl_grid );
free( init_file );
free( egrid_file );
ecl_grid_free( ecl_grid );
}
#endif
// Read an imported fortio data file using Ert.
// Data stored in 'integer_field_map_' and 'floating_field_map_'.
void EclipseGridParser::getNumericErtFields(const string& filename)
{
#ifdef HAVE_ERT
// Read file
ecl_file_type * ecl_file = ecl_file_open(filename.c_str());
if (ecl_file == NULL) {
THROW("Could not open IMPORTed file " << filename);
}
const int num_kw = ecl_file_get_size(ecl_file);
std::vector<double> double_vec;
std::vector<int> int_vec;
for (int i=0; i<num_kw; ++i) {
ecl_kw_type * ecl_kw = ecl_file_iget_kw(ecl_file, i);
const char* keyword = ecl_kw_get_header(ecl_kw);
FieldType field_type = classifyKeyword(keyword);
if (field_type == Unknown) {
ignored_fields_.insert(keyword);
cout << "*** Warning: keyword " << keyword << " is unknown." << endl;
continue;
} else {
#ifdef VERBOSE
cout << "Imported keyword found: " << keyword << endl;
#endif
}
ecl_type_enum ecl_type = ecl_kw_get_type(ecl_kw);
int data_size = ecl_kw_get_size(ecl_kw);
switch(ecl_type) {
case ECL_FLOAT_TYPE : {
double_vec.resize(data_size);
ecl_kw_get_data_as_double(ecl_kw, &double_vec[0]);
floating_field_map_[keyword] = double_vec;
break;
}
case ECL_DOUBLE_TYPE : {
double_vec.resize(data_size);
ecl_kw_get_memcpy_double_data(ecl_kw, &double_vec[0]);
floating_field_map_[keyword] = double_vec;
break;
}
case ECL_INT_TYPE : {
int_vec.resize(data_size);
ecl_kw_get_memcpy_int_data(ecl_kw, &int_vec[0]);
integer_field_map_[keyword] = int_vec;
break;
}
default: {
std::cout << "Ignored non-numeric type in file: " << filename << " Keyword="
<< keyword << " Size=" << ecl_kw_get_size(ecl_kw)
<< " Type=" << ecl_util_get_type_name(ecl_kw_get_type(ecl_kw))
<< std::endl;
break;
}
}
}
ecl_file_close(ecl_file);
#else
THROW("Cannot use IMPORT keyword without ert library support. Reconfigure opm-core with --with-ert and recompile.");
#endif // HAVE_ERT
}
} // namespace Opm } // namespace Opm

24
opm/core/eclipse/EclipseGridParser.hpp
View File

@ -45,6 +45,13 @@ along with OpenRS. If not, see <http://www.gnu.org/licenses/>.
#include <opm/core/eclipse/EclipseUnits.hpp> #include <opm/core/eclipse/EclipseUnits.hpp>
#include <opm/core/utility/Factory.hpp> #include <opm/core/utility/Factory.hpp>
#include <opm/core/grid/cornerpoint_grid.h>
#ifdef HAVE_ERT
#include <ecl_kw.h>
#include <ecl_grid.h>
#endif
namespace Opm namespace Opm
{ {
@ -191,11 +198,28 @@ public:
/// The units specified by the eclipse file read. /// The units specified by the eclipse file read.
const EclipseUnits& units() const; const EclipseUnits& units() const;
struct grdecl get_grdecl() const;
#ifdef HAVE_ERT
void saveEGRID_INIT( const std::string& output_dir , const std::string& basename, bool fmt_file = false);
void saveEGRID( const std::string & filename );
void saveINIT( const std::string & filename , const ecl_grid_type * ecl_grid);
ecl_grid_type * newGrid( );
#endif
private: private:
#ifdef HAVE_ERT
ecl_kw_type * newEclKW(const std::string &keyword , ecl_type_enum ecl_type) const;
void save_kw( fortio_type * fortio , const std::string & kw , ecl_type_enum ecl_type);
#endif
SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname); SpecialFieldPtr createSpecialField(std::istream& is, const std::string& fieldname);
SpecialFieldPtr cloneSpecialField(const std::string& fieldname, SpecialFieldPtr cloneSpecialField(const std::string& fieldname,
const std::tr1::shared_ptr<SpecialBase> original); const std::tr1::shared_ptr<SpecialBase> original);
void readImpl(std::istream& is); void readImpl(std::istream& is);
void getNumericErtFields(const std::string& filename);
std::string directory_; std::string directory_;

3
opm/core/grid/cpgpreprocess/preprocess.c
View File

@ -590,6 +590,9 @@ get_zcorn_sign(int nx, int ny, int nz, const int *actnum,
c1 = i/2 + nx*(j/2 + ny*(k/2)); c1 = i/2 + nx*(j/2 + ny*(k/2));
c2 = i/2 + nx*(j/2 + ny*((k+1)/2)); c2 = i/2 + nx*(j/2 + ny*((k+1)/2));
assert (c1 < (nx * ny * nz));
assert (c2 < (nx * ny * nz));
if (((actnum == NULL) || if (((actnum == NULL) ||
(actnum[c1] && actnum[c2])) (actnum[c1] && actnum[c2]))
&& (z2 < z1)) { && (z2 < z1)) {

1
opm/core/grid/cpgpreprocess/preprocess.h
View File

@ -58,6 +58,7 @@ extern "C" {
const double *coord; /**< Pillar end-points. */ const double *coord; /**< Pillar end-points. */
const double *zcorn; /**< Corner-point depths. */ const double *zcorn; /**< Corner-point depths. */
const int *actnum; /**< Explicit "active" map. May be NULL.*/ const int *actnum; /**< Explicit "active" map. May be NULL.*/
const double *mapaxes; /**< 6 Element rotation vector - can be NULL. */
}; };
/** /**

32
opm/core/utility/DataMap.hpp Normal file
View File

@ -0,0 +1,32 @@
/*
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_DATAMAP_HEADER_INCLUDED
#define OPM_DATAMAP_HEADER_INCLUDED
#include <string>
#include <map>
namespace Opm
{
/// Intended to map strings (giving the output field names) to data.
typedef std::map<std::string, const std::vector<double>*> DataMap;
}
#endif

108
opm/core/utility/writeECLData.cpp Normal file
View File

@ -0,0 +1,108 @@
/*
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 <opm/core/grid.h>
#include <opm/core/simulator/SimulatorTimer.hpp>
#include <opm/core/utility/writeECLData.hpp>
#include <opm/core/utility/Units.hpp>
#include <vector>
#include <ecl_grid.h>
#include <ecl_util.h>
#include <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 ) {
ecl_kw_type * ecl_kw = ecl_kw_alloc( kw_name.c_str() , data->size() / stride , ECL_FLOAT_TYPE );
if (grid.global_cell == NULL) {
for (int i=0; i < grid.number_of_cells; i++)
ecl_kw_iset_float( ecl_kw , i , (*data)[i*stride + offset]);
} else {
for (int i=0; i < grid.number_of_cells; i++)
ecl_kw_iset_float( ecl_kw , grid.global_cell[i] , (*data)[i*stride + offset]);
}
return ecl_kw;
}
void writeECLData(const UnstructuredGrid& grid,
const DataMap& data,
const SimulatorTimer& simtimer,
const std::string& output_dir,
const std::string& base_name) {
int step = simtimer.currentStepNum();
ecl_file_enum file_type = ECL_UNIFIED_RESTART_FILE;
bool fmt_file = true;
char * filename = ecl_util_alloc_filename(output_dir.c_str() , base_name.c_str() , file_type , fmt_file , step );
int phases = ECL_OIL_PHASE + ECL_WATER_PHASE;
double days = Opm::unit::convert::to(simtimer.currentTime(), Opm::unit::day);
time_t date = 0;
int nx = grid.cartdims[0];
int ny = grid.cartdims[1];
int nz = grid.cartdims[2];
int nactive = grid.number_of_cells;
ecl_rst_file_type * rst_file;
if (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_rst_file_fwrite_header( rst_file , step , date , days , nx , ny , nz , nactive , phases );
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()) {
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);
}
}

47
opm/core/utility/writeECLData.hpp Normal file
View File

@ -0,0 +1,47 @@
/*
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 <string>
#include <map>
#include <vector>
#include <tr1/array>
#include <iosfwd>
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
struct UnstructuredGrid;
namespace Opm
{
// ECLIPSE output for general grids.
void writeECLData(const UnstructuredGrid& grid,
const DataMap& data,
const SimulatorTimer& simtimer,
const std::string& output_dir,
const std::string& base_name);
}
#endif

501
opm/core/utility/writeVtkData.cpp
View File

@ -18,6 +18,7 @@
*/ */
#include <opm/core/utility/writeVtkData.hpp> #include <opm/core/utility/writeVtkData.hpp>
#include <opm/core/utility/DataMap.hpp>
#include <opm/core/grid.h> #include <opm/core/grid.h>
#include <opm/core/utility/ErrorMacros.hpp> #include <opm/core/utility/ErrorMacros.hpp>
#include <boost/lexical_cast.hpp> #include <boost/lexical_cast.hpp>
@ -34,63 +35,63 @@ namespace Opm
{ {
void writeVtkData(const std::tr1::array<int, 3>& dims, void writeVtkData(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size, const std::tr1::array<double, 3>& cell_size,
const DataMap& data, const DataMap& data,
std::ostream& os) std::ostream& os)
{ {
// Dimension is hardcoded in the prototype and the next two lines, // Dimension is hardcoded in the prototype and the next two lines,
// but the rest is flexible (allows dimension == 2 or 3). // but the rest is flexible (allows dimension == 2 or 3).
int dimension = 3; int dimension = 3;
int num_cells = dims[0]*dims[1]*dims[2]; int num_cells = dims[0]*dims[1]*dims[2];
ASSERT(dimension == 2 || dimension == 3); ASSERT(dimension == 2 || dimension == 3);
ASSERT(num_cells == dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2])); ASSERT(num_cells == dims[0]*dims[1]* (dimension == 2 ? 1 : dims[2]));
os << "# vtk DataFile Version 2.0\n"; os << "# vtk DataFile Version 2.0\n";
os << "Structured Grid\n \n"; os << "Structured Grid\n \n";
os << "ASCII \n"; os << "ASCII \n";
os << "DATASET STRUCTURED_POINTS\n"; os << "DATASET STRUCTURED_POINTS\n";
os << "DIMENSIONS " os << "DIMENSIONS "
<< dims[0] + 1 << " " << dims[0] + 1 << " "
<< dims[1] + 1 << " "; << dims[1] + 1 << " ";
if (dimension == 3) { if (dimension == 3) {
os << dims[2] + 1; os << dims[2] + 1;
} else { } else {
os << 1; os << 1;
} }
os << "\n"; os << "\n";
os << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n"; os << "ORIGIN " << 0.0 << " " << 0.0 << " " << 0.0 << "\n";
os << "SPACING " << cell_size[0] << " " << cell_size[1]; os << "SPACING " << cell_size[0] << " " << cell_size[1];
if (dimension == 3) { if (dimension == 3) {
os << " " << cell_size[2]; os << " " << cell_size[2];
} else { } else {
os << " " << 0.0; os << " " << 0.0;
} }
os << "\n"; os << "\n";
os << "\nCELL_DATA " << num_cells << '\n'; os << "\nCELL_DATA " << num_cells << '\n';
for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
std::string name = dit->first; std::string name = dit->first;
os << "SCALARS " << name << " float" << '\n'; os << "SCALARS " << name << " float" << '\n';
os << "LOOKUP_TABLE " << name << "_table " << '\n'; os << "LOOKUP_TABLE " << name << "_table " << '\n';
const std::vector<double>& field = *(dit->second); const std::vector<double>& field = *(dit->second);
// We always print only the first data item for every // We always print only the first data item for every
// cell, using 'stride'. // cell, using 'stride'.
// This is a hack to get water saturation nicely. // This is a hack to get water saturation nicely.
// \TODO: Extend to properly printing vector data. // \TODO: Extend to properly printing vector data.
const int stride = field.size()/num_cells; const int stride = field.size()/num_cells;
const int num_per_line = 5; const int num_per_line = 5;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
os << field[stride*c] << ' '; os << field[stride*c] << ' ';
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
} }
typedef std::map<std::string, std::string> PMap; typedef std::map<std::string, std::string> PMap;
@ -98,219 +99,219 @@ namespace Opm
struct Tag struct Tag
{ {
Tag(const std::string& tag, const PMap& props, std::ostream& os) Tag(const std::string& tag, const PMap& props, std::ostream& os)
: name_(tag), os_(os) : name_(tag), os_(os)
{ {
indent(os); indent(os);
os << "<" << tag; os << "<" << tag;
for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) { for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) {
os << " " << it->first << "=\"" << it->second << "\""; os << " " << it->first << "=\"" << it->second << "\"";
} }
os << ">\n"; os << ">\n";
++indent_; ++indent_;
} }
Tag(const std::string& tag, std::ostream& os) Tag(const std::string& tag, std::ostream& os)
: name_(tag), os_(os) : name_(tag), os_(os)
{ {
indent(os); indent(os);
os << "<" << tag << ">\n"; os << "<" << tag << ">\n";
++indent_; ++indent_;
} }
~Tag() ~Tag()
{ {
--indent_; --indent_;
indent(os_); indent(os_);
os_ << "</" << name_ << ">\n"; os_ << "</" << name_ << ">\n";
} }
static void indent(std::ostream& os) static void indent(std::ostream& os)
{ {
for (int i = 0; i < indent_; ++i) { for (int i = 0; i < indent_; ++i) {
os << " "; os << " ";
} }
} }
private: private:
static int indent_; static int indent_;
std::string name_; std::string name_;
std::ostream& os_; std::ostream& os_;
}; };
int Tag::indent_ = 0; int Tag::indent_ = 0;
void writeVtkData(const UnstructuredGrid& grid, void writeVtkData(const UnstructuredGrid& grid,
const DataMap& data, const DataMap& data,
std::ostream& os) std::ostream& os)
{ {
if (grid.dimensions != 3) { if (grid.dimensions != 3) {
THROW("Vtk output for 3d grids only"); THROW("Vtk output for 3d grids only");
} }
os.precision(12); os.precision(12);
os << "<?xml version=\"1.0\"?>\n"; os << "<?xml version=\"1.0\"?>\n";
PMap pm; PMap pm;
pm["type"] = "UnstructuredGrid"; pm["type"] = "UnstructuredGrid";
Tag vtkfiletag("VTKFile", pm, os); Tag vtkfiletag("VTKFile", pm, os);
Tag ugtag("UnstructuredGrid", os); Tag ugtag("UnstructuredGrid", os);
int num_pts = grid.number_of_nodes; int num_pts = grid.number_of_nodes;
int num_cells = grid.number_of_cells; int num_cells = grid.number_of_cells;
pm.clear(); pm.clear();
pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts); pm["NumberOfPoints"] = boost::lexical_cast<std::string>(num_pts);
pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells); pm["NumberOfCells"] = boost::lexical_cast<std::string>(num_cells);
Tag piecetag("Piece", pm, os); Tag piecetag("Piece", pm, os);
{ {
Tag pointstag("Points", os); Tag pointstag("Points", os);
pm.clear(); pm.clear();
pm["type"] = "Float64"; pm["type"] = "Float64";
pm["Name"] = "Coordinates"; pm["Name"] = "Coordinates";
pm["NumberOfComponents"] = "3"; pm["NumberOfComponents"] = "3";
pm["format"] = "ascii"; pm["format"] = "ascii";
Tag datag("DataArray", pm, os); Tag datag("DataArray", pm, os);
for (int i = 0; i < num_pts; ++i) { for (int i = 0; i < num_pts; ++i) {
Tag::indent(os); Tag::indent(os);
os << grid.node_coordinates[3*i + 0] << ' ' os << grid.node_coordinates[3*i + 0] << ' '
<< grid.node_coordinates[3*i + 1] << ' ' << grid.node_coordinates[3*i + 1] << ' '
<< grid.node_coordinates[3*i + 2] << '\n'; << grid.node_coordinates[3*i + 2] << '\n';
} }
} }
{ {
Tag cellstag("Cells", os); Tag cellstag("Cells", os);
pm.clear(); pm.clear();
pm["type"] = "Int32"; pm["type"] = "Int32";
pm["NumberOfComponents"] = "1"; pm["NumberOfComponents"] = "1";
pm["format"] = "ascii"; pm["format"] = "ascii";
std::vector<int> cell_numpts; std::vector<int> cell_numpts;
cell_numpts.reserve(num_cells); cell_numpts.reserve(num_cells);
{ {
pm["Name"] = "connectivity"; pm["Name"] = "connectivity";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
int hf = 0; int hf = 0;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
std::set<int> cell_pts; std::set<int> cell_pts;
for (; hf < grid.cell_facepos[c+1]; ++hf) { for (; hf < grid.cell_facepos[c+1]; ++hf) {
int f = grid.cell_faces[hf]; int f = grid.cell_faces[hf];
const int* fnbeg = grid.face_nodes + grid.face_nodepos[f]; const int* fnbeg = grid.face_nodes + grid.face_nodepos[f];
const int* fnend = grid.face_nodes + grid.face_nodepos[f+1]; const int* fnend = grid.face_nodes + grid.face_nodepos[f+1];
cell_pts.insert(fnbeg, fnend); cell_pts.insert(fnbeg, fnend);
} }
cell_numpts.push_back(cell_pts.size()); cell_numpts.push_back(cell_pts.size());
Tag::indent(os); Tag::indent(os);
std::copy(cell_pts.begin(), cell_pts.end(), std::copy(cell_pts.begin(), cell_pts.end(),
std::ostream_iterator<int>(os, " ")); std::ostream_iterator<int>(os, " "));
os << '\n'; os << '\n';
} }
} }
{ {
pm["Name"] = "offsets"; pm["Name"] = "offsets";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
int offset = 0; int offset = 0;
const int num_per_line = 10; const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) { if (c % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
offset += cell_numpts[c]; offset += cell_numpts[c];
os << offset << ' '; os << offset << ' ';
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
std::vector<int> cell_foffsets; std::vector<int> cell_foffsets;
cell_foffsets.reserve(num_cells); cell_foffsets.reserve(num_cells);
{ {
pm["Name"] = "faces"; pm["Name"] = "faces";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
const int* fp = grid.cell_facepos; const int* fp = grid.cell_facepos;
int offset = 0; int offset = 0;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
Tag::indent(os); Tag::indent(os);
os << fp[c+1] - fp[c] << '\n'; os << fp[c+1] - fp[c] << '\n';
++offset; ++offset;
for (int hf = fp[c]; hf < fp[c+1]; ++hf) { for (int hf = fp[c]; hf < fp[c+1]; ++hf) {
int f = grid.cell_faces[hf]; int f = grid.cell_faces[hf];
const int* np = grid.face_nodepos; const int* np = grid.face_nodepos;
int f_num_pts = np[f+1] - np[f]; int f_num_pts = np[f+1] - np[f];
Tag::indent(os); Tag::indent(os);
os << f_num_pts << ' '; os << f_num_pts << ' ';
++offset; ++offset;
std::copy(grid.face_nodes + np[f], std::copy(grid.face_nodes + np[f],
grid.face_nodes + np[f+1], grid.face_nodes + np[f+1],
std::ostream_iterator<int>(os, " ")); std::ostream_iterator<int>(os, " "));
os << '\n'; os << '\n';
offset += f_num_pts; offset += f_num_pts;
} }
cell_foffsets.push_back(offset); cell_foffsets.push_back(offset);
} }
} }
{ {
pm["Name"] = "faceoffsets"; pm["Name"] = "faceoffsets";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
const int num_per_line = 10; const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) { if (c % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
os << cell_foffsets[c] << ' '; os << cell_foffsets[c] << ' ';
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
{ {
pm["type"] = "UInt8"; pm["type"] = "UInt8";
pm["Name"] = "types"; pm["Name"] = "types";
Tag t("DataArray", pm, os); Tag t("DataArray", pm, os);
const int num_per_line = 10; const int num_per_line = 10;
for (int c = 0; c < num_cells; ++c) { for (int c = 0; c < num_cells; ++c) {
if (c % num_per_line == 0) { if (c % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
os << "42 "; os << "42 ";
if (c % num_per_line == num_per_line - 1 if (c % num_per_line == num_per_line - 1
|| c == num_cells - 1) { || c == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
} }
{ {
pm.clear(); pm.clear();
if (data.find("saturation") != data.end()) { if (data.find("saturation") != data.end()) {
pm["Scalars"] = "saturation"; pm["Scalars"] = "saturation";
} else if (data.find("pressure") != data.end()) { } else if (data.find("pressure") != data.end()) {
pm["Scalars"] = "pressure"; pm["Scalars"] = "pressure";
} }
Tag celldatatag("CellData", pm, os); Tag celldatatag("CellData", pm, os);
pm.clear(); pm.clear();
pm["NumberOfComponents"] = "1"; pm["NumberOfComponents"] = "1";
pm["format"] = "ascii"; pm["format"] = "ascii";
pm["type"] = "Float64"; pm["type"] = "Float64";
for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) {
pm["Name"] = dit->first; pm["Name"] = dit->first;
const std::vector<double>& field = *(dit->second); const std::vector<double>& field = *(dit->second);
const int num_comps = field.size()/grid.number_of_cells; const int num_comps = field.size()/grid.number_of_cells;
pm["NumberOfComponents"] = boost::lexical_cast<std::string>(num_comps); pm["NumberOfComponents"] = boost::lexical_cast<std::string>(num_comps);
Tag ptag("DataArray", pm, os); Tag ptag("DataArray", pm, os);
const int num_per_line = num_comps == 1 ? 5 : num_comps; const int num_per_line = num_comps == 1 ? 5 : num_comps;
for (int item = 0; item < num_cells*num_comps; ++item) { for (int item = 0; item < num_cells*num_comps; ++item) {
if (item % num_per_line == 0) { if (item % num_per_line == 0) {
Tag::indent(os); Tag::indent(os);
} }
double value = field[item]; double value = field[item];
if (std::fabs(value) < std::numeric_limits<double>::min()) { if (std::fabs(value) < std::numeric_limits<double>::min()) {
// Avoiding denormal numbers to work around // Avoiding denormal numbers to work around
// bug in Paraview. // bug in Paraview.
value = 0.0; value = 0.0;
} }
os << value << ' '; os << value << ' ';
if (item % num_per_line == num_per_line - 1 if (item % num_per_line == num_per_line - 1
|| item == num_cells - 1) { || item == num_cells - 1) {
os << '\n'; os << '\n';
} }
} }
} }
} }
} }
} // namespace Opm } // namespace Opm

14
opm/core/utility/writeVtkData.hpp
View File

@ -26,25 +26,23 @@
#include <vector> #include <vector>
#include <tr1/array> #include <tr1/array>
#include <iosfwd> #include <iosfwd>
#include <opm/core/utility/DataMap.hpp>
struct UnstructuredGrid; struct UnstructuredGrid;
namespace Opm namespace Opm
{ {
/// Intended to map strings (giving the output field names) to data.
typedef std::map<std::string, const std::vector<double>*> DataMap;
/// Vtk output for cartesian grids. /// Vtk output for cartesian grids.
void writeVtkData(const std::tr1::array<int, 3>& dims, void writeVtkData(const std::tr1::array<int, 3>& dims,
const std::tr1::array<double, 3>& cell_size, const std::tr1::array<double, 3>& cell_size,
const DataMap& data, const DataMap& data,
std::ostream& os); std::ostream& os);
/// Vtk output for general grids. /// Vtk output for general grids.
void writeVtkData(const UnstructuredGrid& grid, void writeVtkData(const UnstructuredGrid& grid,
const DataMap& data, const DataMap& data,
std::ostream& os); std::ostream& os);
} // namespace Opm } // namespace Opm
#endif // OPM_WRITEVTKDATA_HEADER_INCLUDED #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED

11
tests/Makefile.am
View File

@ -1,8 +1,9 @@
AM_CPPFLAGS = \ AM_CPPFLAGS = \
-I$(top_srcdir) \ -I$(top_srcdir) \
$(OPM_BOOST_CPPFLAGS) ${SUPERLU_CPPFLAGS} $(OPM_BOOST_CPPFLAGS) ${SUPERLU_CPPFLAGS} \
$(ERT_CPPFLAGS)
AM_LDFLAGS = $(OPM_BOOST_LDFLAGS) AM_LDFLAGS = $(OPM_BOOST_LDFLAGS) $(ERT_LDFLAGS)
LDADD = $(top_builddir)/lib/libopmcore.la LDADD = $(top_builddir)/lib/libopmcore.la
@ -77,3 +78,9 @@ if BUILD_AGMG
noinst_PROGRAMS += test_agmg noinst_PROGRAMS += test_agmg
test_agmg_SOURCES = test_agmg.cpp test_agmg_SOURCES = test_agmg.cpp
endif endif
if HAVE_ERT
noinst_PROGRAMS += test_ert
test_ert_SOURCES = test_ert.cpp
endif

260
tests/test_ert.cpp Normal file
View File

@ -0,0 +1,260 @@
/*
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/eclipse/EclipseGridParser.hpp>
#include <opm/core/GridManager.hpp>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid.h>
#include <cstdio>
#include <boost/scoped_ptr.hpp>
#include <opm/core/fluid/IncompPropertiesBasic.hpp>
#include <opm/core/fluid/IncompPropertiesFromDeck.hpp>
#include <ecl_grid.h>
using namespace std;
#if 0
static
void cell_nodes(const UnstructuredGrid * c_grid , int cell , std::vector<int>& nodes) {
int face_offset = c_grid->cell_facepos[cell];
int num_faces = c_grid->cell_facepos[cell + 1] - face_offset;
nodes.clear();
//printf("cell: %d \n",cell);
for (int iface = 0; iface < num_faces; iface++) {
int face = c_grid->cell_faces[ face_offset + iface];
//printf("face[%d] = %d \n",iface , face );
{
int node_offset = c_grid->face_nodepos[ face ];
int num_nodes = c_grid->face_nodepos[ face + 1] - node_offset;
for (int inode = 0; inode < num_nodes; inode++) {
int node = c_grid->face_nodes[ inode + node_offset ];
//printf(" node[%d] = %d \n",inode , node);
nodes.push_back( node );
}
}
}
{
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
sort( nodes.begin() , nodes.end());
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
unique( nodes.begin() , nodes.end() );
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
nodes.resize( 8 );
/*for (int i =0; i < nodes.size(); i++)
std::cout << nodes[i] << " ";
std::cout << "\n";
*/
}
}
#endif
/*
static
void eclExport(Opm::GridManager& grid) {
const UnstructuredGrid * c_grid = grid.c_grid();
printf("dimensions : %d \n",c_grid->dimensions);
printf("number of cells : %d \n",c_grid->number_of_cells);
printf("number of nodes : %d \n",c_grid->number_of_nodes);
printf("number of faces : %d \n",c_grid->number_of_faces);
printf("length(face_nodes) : %d \n",c_grid->face_nodepos[ c_grid->number_of_faces ]);
printf("cartdims : %d %d %d \n",
c_grid->cartdims[0] ,
c_grid->cartdims[1] ,
c_grid->cartdims[2]);
printf("global_cell : %d %d %d %d %d\n",
c_grid->global_cell[0] ,
c_grid->global_cell[1] ,
c_grid->global_cell[2] ,
c_grid->global_cell[3] ,
c_grid->global_cell[4]);
{
std::vector<int> nodes;
cell_nodes( c_grid , 10 , nodes );
cell_nodes( c_grid , 15 , nodes );
cell_nodes( c_grid , 20 , nodes );
cell_nodes( c_grid , 25 , nodes );
}
{
ecl_grid_type * ecl_grid;
int num_coords = c_grid->number_of_cells;
int coords_size = 4;
int nx = c_grid->cartdims[0];
int ny = c_grid->cartdims[1];
int nz = c_grid->cartdims[2];
int ** coords;
float ** corners;
// float * mapaxes = NULL;
std::vector<int> nodes;
corners = (float **) malloc( num_coords * sizeof * corners );
coords = (int **) malloc( num_coords * sizeof * coords );
{
int c;
for (c=0; c < num_coords; c++) {
corners[c] = (float *) malloc( 24 * sizeof * corners[c] );
coords[c] = (int *) malloc( coords_size * sizeof * coords[c] );
}
for (c=0; c < num_coords; c++) {
cell_nodes( c_grid , c , nodes );
for (int p=0; p < 8; p++) {
int n = nodes[p];
for (int d=0; d < c_grid->dimensions; d++)
corners[c][3*p + d] = c_grid->node_coordinates[ c_grid->dimensions * n + d ];
}
{
int i,j,k;
{
int g = c_grid->global_cell[ c ];
k = g / nx*ny; g -= k * nx*ny;
j = g / nx; g -= j * nx;
i = g;
}
coords[c][0] = i + 1;
coords[c][1] = j + 1;
coords[c][2] = k + 1;
coords[c][3] = c_grid->global_cell[ c ] + 1;
}
}
}
ecl_grid = ecl_grid_alloc( "/private/joaho/ERT/NR/libenkf/src/Gurbat/EXAMPLE_01_BASE.EGRID" );
printf("Grid loaded ... \n");
ecl_grid_free( ecl_grid );
printf("Grid discarded ... \n");
ecl_grid = ecl_grid_alloc_GRID_data( num_coords , nx , ny , nz , coords_size , coords , corners , NULL );
ecl_grid_fwrite_GRID( ecl_grid , "/tmp/test.GRID" );
{
FILE * stream = fopen( "/tmp/test.grdecl" , "w");
ecl_grid_fprintf_grdecl( ecl_grid , stream );
fclose( stream );
}
ecl_grid_free( ecl_grid );
{
for (int c=0; c < num_coords; c++) {
free(corners[c]);
free(coords[c]);
}
}
free( corners );
free( coords );
}
}
*/
/*
#ifdef HAVE_ERT
ecl_grid_type * create_ecl_grid( const struct UnstructuredGrid * g) {
int num_coords = g->number_of_cells;
int nx = g->cartdims[0];
int ny = g->cartdims[1];
int nz = g->cartdims[2];
int coords_size = 4;
int ** coords;
float ** corners;
float * mapaxes = NULL;
corners = malloc( num_coords * sizeof * corners );
coords = malloc( num_coords * sizeof * coords );
{
for (int c=0; c < num_coords; c++) {
corners[c] = malloc( 24 * sizeof * corners[0] );
coords[c] = malloc( coords_size * sizeof * coords[0] );
}
}
{
for (int k=0; k < nz; k++) {
for (int j=0; j < ny; j++) {
for (int i=0; i < nx; i++) {
int global_index = i + j*nx + k*nx*ny;
coords[global_index][0] = i;
coords[global_index][1] = j;
coords[global_index][2] = k;
coords[global_index][3] = 1;
}
}
}
}
{
for (int c=0; c < num_coords; c++) {
free(corners[c]);
free(coords[c]);
}
}
free( corners );
free( coords );
}
#endif
*/
// struct grdecl : opm/core/grid/cpgpreprocess/preprocess.h
// struct processes_grid : opm/core/grid/cpgpreprocess/preprocess.h
int main(int /*argc*/ , char **argv)
{
std::string filename( argv[1] );
boost::scoped_ptr<Opm::GridManager> grid;
boost::scoped_ptr<Opm::IncompPropertiesInterface> props;
Opm::EclipseGridParser eclParser(filename , false);
//eclParser.saveEGRID_INIT("/tmp" , "OPM" );
grid.reset(new Opm::GridManager(eclParser));
props.reset(new Opm::IncompPropertiesFromDeck(eclParser , *grid->c_grid()));
}