diff --git a/Makefile.am b/Makefile.am index 70a16dfd..787af5f6 100644 --- a/Makefile.am +++ b/Makefile.am @@ -12,6 +12,7 @@ lib_LTLIBRARIES = lib/libopmcore.la # Build-time flags needed to build libopmcore.la AM_CPPFLAGS = \ +$(ERT_CPPFLAGS) \ $(OPM_BOOST_CPPFLAGS) \ ${SUPERLU_CPPFLAGS} @@ -22,6 +23,7 @@ ${SUPERLU_CPPFLAGS} lib_libopmcore_la_LDFLAGS = \ -R $(OPM_BOOST_LIBDIR) \ $(OPM_BOOST_LDFLAGS) \ +$(ERT_LDFLAGS) ${SUPERLU_LDFLAGS} lib_libopmcore_la_LIBADD = \ @@ -29,6 +31,7 @@ $(BOOST_FILESYSTEM_LIB) \ $(BOOST_SYSTEM_LIB) \ $(BOOST_DATE_TIME_LIB) \ $(BOOST_UNIT_TEST_FRAMEWORK_LIB) \ +$(ERT_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/utility/Average.hpp \ opm/core/utility/ColumnExtract.hpp \ +opm/core/utility/DataMap.hpp \ opm/core/utility/ErrorMacros.hpp \ opm/core/utility/Factory.hpp \ opm/core/utility/MonotCubicInterpolator.hpp \ @@ -292,6 +296,15 @@ opm/core/linalg/LinearSolverUmfpack.hpp 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 lib_libopmcore_la_SOURCES += \ opm/core/linalg/LinearSolverIstl.cpp diff --git a/configure.ac b/configure.ac index 045faeda..4644d349 100644 --- a/configure.ac +++ b/configure.ac @@ -39,6 +39,8 @@ OPM_CORE_CHECKS OPM_DYNLINK_BOOST_TEST +ERT + dnl Substitute Autoconf's abs_*dir variables into the Makefiles for the dnl benefit of external code that uses these variables to derive dnl locations (e.g., Dune's DUNE_CHECK_MODULES macro). Automakes prior diff --git a/examples/Makefile.am b/examples/Makefile.am index 1f068bcf..145fbbdd 100644 --- a/examples/Makefile.am +++ b/examples/Makefile.am @@ -1,6 +1,9 @@ # Build-time flags needed to form example programs +ERT_INCLUDE_PATH = $(ERT_ROOT)/include + AM_CPPFLAGS = \ -I$(top_srcdir) \ +-I$(ERT_INCLUDE_PATH) \ $(OPM_BOOST_CPPFLAGS) # All targets link to the library diff --git a/examples/spu_2p.cpp b/examples/spu_2p.cpp index 25cd3f6f..755ef202 100644 --- a/examples/spu_2p.cpp +++ b/examples/spu_2p.cpp @@ -94,14 +94,19 @@ #include +#ifdef HAVE_ERT +#include +#endif + static void outputState(const UnstructuredGrid& grid, const Opm::TwophaseState& state, - const int step, + const Opm::SimulatorTimer& simtimer, const std::string& output_dir) { // Write data in VTK format. + int step = simtimer.currentStepNum(); std::ostringstream vtkfilename; vtkfilename << output_dir << "/output-" << std::setw(3) << std::setfill('0') << step << ".vtu"; std::ofstream vtkfile(vtkfilename.str().c_str()); @@ -115,6 +120,9 @@ static void outputState(const UnstructuredGrid& grid, Opm::estimateCellVelocity(grid, state.faceflux(), cell_velocity); dm["velocity"] = &cell_velocity; Opm::writeVtkData(grid, dm, vtkfile); +#ifdef HAVE_ERT + Opm::writeECLData(grid , dm , simtimer , output_dir , "OPM" ); +#endif // Write data (not grid) in Matlab format 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. simtimer.report(std::cout); 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. @@ -699,7 +707,7 @@ main(int argc, char** argv) << "\n Transport time: " << ttime << std::endl; if (output) { - outputState(*grid->c_grid(), state, simtimer.currentStepNum(), output_dir); + outputState(*grid->c_grid(), state, simtimer, output_dir); outputWaterCut(watercut, output_dir); if (wells->c_wells()) { outputWellReport(wellreport, output_dir); diff --git a/m4/ert.m4 b/m4/ert.m4 new file mode 100644 index 00000000..28067892 --- /dev/null +++ b/m4/ert.m4 @@ -0,0 +1,67 @@ +AC_DEFUN([_ERT_SOURCE_TEXT], +[ +AC_LANG_PROGRAM( +[[ +#include +#include +]],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=], [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} **** +# ]) +]) diff --git a/opm/core/GridManager.cpp b/opm/core/GridManager.cpp index 04d8cd04..d286e504 100644 --- a/opm/core/GridManager.cpp +++ b/opm/core/GridManager.cpp @@ -37,7 +37,8 @@ namespace Opm { // We accept two different ways to specify the grid. // 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, // DEPTHZ and TOPS are not present. // 2. Tensor grid format. @@ -119,29 +120,8 @@ namespace Opm void GridManager::initFromDeckCornerpoint(const Opm::EclipseGridParser& deck) { // Extract data from deck. - const std::vector& zcorn = deck.getFloatingPointValue("ZCORN"); - const std::vector& coord = deck.getFloatingPointValue("COORD"); - const int* actnum = 0; - if (deck.hasField("ACTNUM")) { - actnum = &(deck.getIntegerValue("ACTNUM")[0]); - } - std::vector 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. - struct grdecl 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]; + struct grdecl grdecl = deck.get_grdecl(); // Process grid. ug_ = create_grid_cornerpoint(&grdecl, 0.0); diff --git a/opm/core/eclipse/EclipseGridParser.cpp b/opm/core/eclipse/EclipseGridParser.cpp index b4ea9ee3..746794db 100644 --- a/opm/core/eclipse/EclipseGridParser.cpp +++ b/opm/core/eclipse/EclipseGridParser.cpp @@ -47,8 +47,19 @@ #include #include #include -#include #include +#include +#include + +#ifdef HAVE_ERT +#include +#include +#include +#include +#include +#include +#include +#endif using namespace std; @@ -86,7 +97,7 @@ namespace EclipseKeywords string("MULTPV"), string("PRESSURE"), string("SGAS"), string("SWAT"), string("SOIL"), string("RS"), 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]); @@ -114,7 +125,7 @@ namespace EclipseKeywords const int num_special_fields = sizeof(special_fields) / sizeof(special_fields[0]); string ignore_with_data[] = - { string("MAPUNITS"), string("MAPAXES"), string("GRIDUNIT"), + { string("MAPUNITS"), string("GRIDUNIT"), string("NTG"), string("REGDIMS"), string("WELLDIMS"), string("NSTACK"), string("SATNUM"), string("RPTRST"), string("ROIP"), string("RWIP"), @@ -141,6 +152,9 @@ namespace EclipseKeywords string include_keywords[] = { string("INCLUDE") }; 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 @@ -154,6 +168,7 @@ namespace { IgnoreWithData, IgnoreNoData, Include, + Import, Unknown }; @@ -174,6 +189,8 @@ namespace { return IgnoreNoData; } else if (count(include_keywords, include_keywords + num_include_keywords, keyword)) { return Include; + } else if (count(import_keywords, import_keywords + num_import_keywords, keyword)) { + return Import; } else { return Unknown; } @@ -488,6 +505,14 @@ void EclipseGridParser::readImpl(istream& is) // is >> ignoreSlashLine; break; } + case Import: { + string import_filename = readString(is); + if (!directory_.empty()) { + import_filename = directory_ + '/' + import_filename; + } + getNumericErtFields(import_filename); + break; + } case Unknown: default: ignored_fields_.insert(keyword); @@ -542,6 +567,9 @@ void EclipseGridParser::convertToSI() do_convert = false; // Dimensionless keywords... } else if (key == "PRESSURE") { unit = units_.pressure; + } else if (key == "MAPAXES") { + MESSAGE("Not applying units to MAPAXES yet!"); + unit = 1.0; } else { THROW("Units for field " << key << " not specified. Cannon convert to SI."); } @@ -800,5 +828,310 @@ void EclipseGridParser::computeUnits() THROW("Unknown unit family " << unit_family); } } + + +struct grdecl EclipseGridParser::get_grdecl() const { + struct grdecl grdecl; + + // Extract data from deck. + const std::vector& zcorn = getFloatingPointValue("ZCORN"); + const std::vector& coord = getFloatingPointValue("COORD"); + const int* actnum = NULL; + if (hasField("ACTNUM")) { + actnum = &(getIntegerValue("ACTNUM")[0]); + } + + std::vector 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 &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 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 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::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_vec; + std::vector int_vec; + for (int i=0; i. #include #include +#include +#ifdef HAVE_ERT +#include +#include +#endif + + namespace Opm { @@ -191,11 +198,28 @@ public: /// The units specified by the eclipse file read. 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: + +#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 cloneSpecialField(const std::string& fieldname, const std::tr1::shared_ptr original); void readImpl(std::istream& is); + void getNumericErtFields(const std::string& filename); std::string directory_; diff --git a/opm/core/grid/cpgpreprocess/preprocess.c b/opm/core/grid/cpgpreprocess/preprocess.c index 84e9e8fe..dc274160 100644 --- a/opm/core/grid/cpgpreprocess/preprocess.c +++ b/opm/core/grid/cpgpreprocess/preprocess.c @@ -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)); c2 = i/2 + nx*(j/2 + ny*((k+1)/2)); + assert (c1 < (nx * ny * nz)); + assert (c2 < (nx * ny * nz)); + if (((actnum == NULL) || (actnum[c1] && actnum[c2])) && (z2 < z1)) { diff --git a/opm/core/grid/cpgpreprocess/preprocess.h b/opm/core/grid/cpgpreprocess/preprocess.h index 13d682bb..a0db38ce 100644 --- a/opm/core/grid/cpgpreprocess/preprocess.h +++ b/opm/core/grid/cpgpreprocess/preprocess.h @@ -58,6 +58,7 @@ extern "C" { const double *coord; /**< Pillar end-points. */ const double *zcorn; /**< Corner-point depths. */ const int *actnum; /**< Explicit "active" map. May be NULL.*/ + const double *mapaxes; /**< 6 Element rotation vector - can be NULL. */ }; /** diff --git a/opm/core/utility/DataMap.hpp b/opm/core/utility/DataMap.hpp new file mode 100644 index 00000000..833f7bcc --- /dev/null +++ b/opm/core/utility/DataMap.hpp @@ -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 . +*/ + +#ifndef OPM_DATAMAP_HEADER_INCLUDED +#define OPM_DATAMAP_HEADER_INCLUDED + +#include +#include + +namespace Opm +{ + /// Intended to map strings (giving the output field names) to data. + typedef std::map*> DataMap; +} + +#endif diff --git a/opm/core/utility/writeECLData.cpp b/opm/core/utility/writeECLData.cpp new file mode 100644 index 00000000..5e9d4205 --- /dev/null +++ b/opm/core/utility/writeECLData.cpp @@ -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 . +*/ + + + +#include +#include +#include +#include + +#include + +#include +#include +#include + + + +namespace Opm +{ + + static ecl_kw_type * ecl_kw_wrapper( const UnstructuredGrid& grid, + const std::string& kw_name , + const std::vector * data , + int offset , + int stride ) { + + 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); + } +} + diff --git a/opm/core/utility/writeECLData.hpp b/opm/core/utility/writeECLData.hpp new file mode 100644 index 00000000..1c963c0d --- /dev/null +++ b/opm/core/utility/writeECLData.hpp @@ -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 . +*/ + +#ifndef OPM_WRITEECLDATA_HEADER_INCLUDED +#define OPM_WRITEECLDATA_HEADER_INCLUDED + + +#include +#include +#include +#include +#include +#include +#include + +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 diff --git a/opm/core/utility/writeVtkData.cpp b/opm/core/utility/writeVtkData.cpp index 07e57172..3c70016e 100644 --- a/opm/core/utility/writeVtkData.cpp +++ b/opm/core/utility/writeVtkData.cpp @@ -18,6 +18,7 @@ */ #include +#include #include #include #include @@ -34,63 +35,63 @@ namespace Opm { void writeVtkData(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const DataMap& data, - std::ostream& os) + const std::tr1::array& cell_size, + const DataMap& data, + std::ostream& os) { - // Dimension is hardcoded in the prototype and the next two lines, - // but the rest is flexible (allows dimension == 2 or 3). - int dimension = 3; - int num_cells = dims[0]*dims[1]*dims[2]; + // 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])); + 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 << "# 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 << "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 << "SPACING " << cell_size[0] << " " << cell_size[1]; + if (dimension == 3) { + os << " " << cell_size[2]; + } else { + os << " " << 0.0; + } + os << "\n"; - os << "\nCELL_DATA " << num_cells << '\n'; - for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { - std::string name = dit->first; - os << "SCALARS " << name << " float" << '\n'; - os << "LOOKUP_TABLE " << name << "_table " << '\n'; - const std::vector& field = *(dit->second); - // We always print only the first data item for every - // cell, using 'stride'. - // This is a hack to get water saturation nicely. - // \TODO: Extend to properly printing vector data. - const int stride = field.size()/num_cells; - const int num_per_line = 5; - for (int c = 0; c < num_cells; ++c) { - os << field[stride*c] << ' '; - if (c % num_per_line == num_per_line - 1 - || c == num_cells - 1) { - os << '\n'; - } - } - } + os << "\nCELL_DATA " << num_cells << '\n'; + for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { + std::string name = dit->first; + os << "SCALARS " << name << " float" << '\n'; + os << "LOOKUP_TABLE " << name << "_table " << '\n'; + const std::vector& field = *(dit->second); + // We always print only the first data item for every + // cell, using 'stride'. + // This is a hack to get water saturation nicely. + // \TODO: Extend to properly printing vector data. + const int stride = field.size()/num_cells; + const int num_per_line = 5; + for (int c = 0; c < num_cells; ++c) { + os << field[stride*c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } } typedef std::map PMap; @@ -98,219 +99,219 @@ namespace Opm struct Tag { - Tag(const std::string& tag, const PMap& props, std::ostream& os) - : name_(tag), os_(os) - { - indent(os); - os << "<" << tag; - for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) { - os << " " << it->first << "=\"" << it->second << "\""; - } - os << ">\n"; - ++indent_; - } - Tag(const std::string& tag, std::ostream& os) - : name_(tag), os_(os) - { - indent(os); - os << "<" << tag << ">\n"; - ++indent_; - } - ~Tag() - { - --indent_; - indent(os_); - os_ << "\n"; - } - static void indent(std::ostream& os) - { - for (int i = 0; i < indent_; ++i) { - os << " "; - } - } + Tag(const std::string& tag, const PMap& props, std::ostream& os) + : name_(tag), os_(os) + { + indent(os); + os << "<" << tag; + for (PMap::const_iterator it = props.begin(); it != props.end(); ++it) { + os << " " << it->first << "=\"" << it->second << "\""; + } + os << ">\n"; + ++indent_; + } + Tag(const std::string& tag, std::ostream& os) + : name_(tag), os_(os) + { + indent(os); + os << "<" << tag << ">\n"; + ++indent_; + } + ~Tag() + { + --indent_; + indent(os_); + os_ << "\n"; + } + static void indent(std::ostream& os) + { + for (int i = 0; i < indent_; ++i) { + os << " "; + } + } private: - static int indent_; - std::string name_; - std::ostream& os_; + static int indent_; + std::string name_; + std::ostream& os_; }; int Tag::indent_ = 0; void writeVtkData(const UnstructuredGrid& grid, - const DataMap& data, - std::ostream& os) + const DataMap& data, + std::ostream& os) { - if (grid.dimensions != 3) { - THROW("Vtk output for 3d grids only"); - } - os.precision(12); - os << "\n"; - PMap pm; - pm["type"] = "UnstructuredGrid"; - Tag vtkfiletag("VTKFile", pm, os); - Tag ugtag("UnstructuredGrid", os); - int num_pts = grid.number_of_nodes; - int num_cells = grid.number_of_cells; - pm.clear(); - pm["NumberOfPoints"] = boost::lexical_cast(num_pts); - pm["NumberOfCells"] = boost::lexical_cast(num_cells); - Tag piecetag("Piece", pm, os); - { - Tag pointstag("Points", os); - pm.clear(); - pm["type"] = "Float64"; - pm["Name"] = "Coordinates"; - pm["NumberOfComponents"] = "3"; - pm["format"] = "ascii"; - Tag datag("DataArray", pm, os); - for (int i = 0; i < num_pts; ++i) { - Tag::indent(os); - os << grid.node_coordinates[3*i + 0] << ' ' - << grid.node_coordinates[3*i + 1] << ' ' - << grid.node_coordinates[3*i + 2] << '\n'; - } - } - { - Tag cellstag("Cells", os); - pm.clear(); - pm["type"] = "Int32"; - pm["NumberOfComponents"] = "1"; - pm["format"] = "ascii"; - std::vector cell_numpts; - cell_numpts.reserve(num_cells); - { - pm["Name"] = "connectivity"; - Tag t("DataArray", pm, os); - int hf = 0; - for (int c = 0; c < num_cells; ++c) { - std::set cell_pts; - for (; hf < grid.cell_facepos[c+1]; ++hf) { - int f = grid.cell_faces[hf]; - const int* fnbeg = grid.face_nodes + grid.face_nodepos[f]; - const int* fnend = grid.face_nodes + grid.face_nodepos[f+1]; - cell_pts.insert(fnbeg, fnend); - } - cell_numpts.push_back(cell_pts.size()); - Tag::indent(os); - std::copy(cell_pts.begin(), cell_pts.end(), - std::ostream_iterator(os, " ")); - os << '\n'; - } - } - { - pm["Name"] = "offsets"; - Tag t("DataArray", pm, os); - int offset = 0; - const int num_per_line = 10; - for (int c = 0; c < num_cells; ++c) { - if (c % num_per_line == 0) { - Tag::indent(os); - } - offset += cell_numpts[c]; - os << offset << ' '; - if (c % num_per_line == num_per_line - 1 - || c == num_cells - 1) { - os << '\n'; - } - } - } - std::vector cell_foffsets; - cell_foffsets.reserve(num_cells); - { - pm["Name"] = "faces"; - Tag t("DataArray", pm, os); - const int* fp = grid.cell_facepos; - int offset = 0; - for (int c = 0; c < num_cells; ++c) { - Tag::indent(os); - os << fp[c+1] - fp[c] << '\n'; - ++offset; - for (int hf = fp[c]; hf < fp[c+1]; ++hf) { - int f = grid.cell_faces[hf]; - const int* np = grid.face_nodepos; - int f_num_pts = np[f+1] - np[f]; - Tag::indent(os); - os << f_num_pts << ' '; - ++offset; - std::copy(grid.face_nodes + np[f], - grid.face_nodes + np[f+1], - std::ostream_iterator(os, " ")); - os << '\n'; - offset += f_num_pts; - } - cell_foffsets.push_back(offset); - } - } - { - pm["Name"] = "faceoffsets"; - Tag t("DataArray", pm, os); - const int num_per_line = 10; - for (int c = 0; c < num_cells; ++c) { - if (c % num_per_line == 0) { - Tag::indent(os); - } - os << cell_foffsets[c] << ' '; - if (c % num_per_line == num_per_line - 1 - || c == num_cells - 1) { - os << '\n'; - } - } - } - { - pm["type"] = "UInt8"; - pm["Name"] = "types"; - Tag t("DataArray", pm, os); - const int num_per_line = 10; - for (int c = 0; c < num_cells; ++c) { - if (c % num_per_line == 0) { - Tag::indent(os); - } - os << "42 "; - if (c % num_per_line == num_per_line - 1 - || c == num_cells - 1) { - os << '\n'; - } - } - } - } - { - pm.clear(); - if (data.find("saturation") != data.end()) { - pm["Scalars"] = "saturation"; - } else if (data.find("pressure") != data.end()) { - pm["Scalars"] = "pressure"; - } - Tag celldatatag("CellData", pm, os); - pm.clear(); - pm["NumberOfComponents"] = "1"; - pm["format"] = "ascii"; - pm["type"] = "Float64"; - for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { - pm["Name"] = dit->first; - const std::vector& field = *(dit->second); - const int num_comps = field.size()/grid.number_of_cells; - pm["NumberOfComponents"] = boost::lexical_cast(num_comps); - Tag ptag("DataArray", pm, os); - const int num_per_line = num_comps == 1 ? 5 : num_comps; - for (int item = 0; item < num_cells*num_comps; ++item) { - if (item % num_per_line == 0) { - Tag::indent(os); - } + if (grid.dimensions != 3) { + THROW("Vtk output for 3d grids only"); + } + os.precision(12); + os << "\n"; + PMap pm; + pm["type"] = "UnstructuredGrid"; + Tag vtkfiletag("VTKFile", pm, os); + Tag ugtag("UnstructuredGrid", os); + int num_pts = grid.number_of_nodes; + int num_cells = grid.number_of_cells; + pm.clear(); + pm["NumberOfPoints"] = boost::lexical_cast(num_pts); + pm["NumberOfCells"] = boost::lexical_cast(num_cells); + Tag piecetag("Piece", pm, os); + { + Tag pointstag("Points", os); + pm.clear(); + pm["type"] = "Float64"; + pm["Name"] = "Coordinates"; + pm["NumberOfComponents"] = "3"; + pm["format"] = "ascii"; + Tag datag("DataArray", pm, os); + for (int i = 0; i < num_pts; ++i) { + Tag::indent(os); + os << grid.node_coordinates[3*i + 0] << ' ' + << grid.node_coordinates[3*i + 1] << ' ' + << grid.node_coordinates[3*i + 2] << '\n'; + } + } + { + Tag cellstag("Cells", os); + pm.clear(); + pm["type"] = "Int32"; + pm["NumberOfComponents"] = "1"; + pm["format"] = "ascii"; + std::vector cell_numpts; + cell_numpts.reserve(num_cells); + { + pm["Name"] = "connectivity"; + Tag t("DataArray", pm, os); + int hf = 0; + for (int c = 0; c < num_cells; ++c) { + std::set cell_pts; + for (; hf < grid.cell_facepos[c+1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* fnbeg = grid.face_nodes + grid.face_nodepos[f]; + const int* fnend = grid.face_nodes + grid.face_nodepos[f+1]; + cell_pts.insert(fnbeg, fnend); + } + cell_numpts.push_back(cell_pts.size()); + Tag::indent(os); + std::copy(cell_pts.begin(), cell_pts.end(), + std::ostream_iterator(os, " ")); + os << '\n'; + } + } + { + pm["Name"] = "offsets"; + Tag t("DataArray", pm, os); + int offset = 0; + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + offset += cell_numpts[c]; + os << offset << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + std::vector cell_foffsets; + cell_foffsets.reserve(num_cells); + { + pm["Name"] = "faces"; + Tag t("DataArray", pm, os); + const int* fp = grid.cell_facepos; + int offset = 0; + for (int c = 0; c < num_cells; ++c) { + Tag::indent(os); + os << fp[c+1] - fp[c] << '\n'; + ++offset; + for (int hf = fp[c]; hf < fp[c+1]; ++hf) { + int f = grid.cell_faces[hf]; + const int* np = grid.face_nodepos; + int f_num_pts = np[f+1] - np[f]; + Tag::indent(os); + os << f_num_pts << ' '; + ++offset; + std::copy(grid.face_nodes + np[f], + grid.face_nodes + np[f+1], + std::ostream_iterator(os, " ")); + os << '\n'; + offset += f_num_pts; + } + cell_foffsets.push_back(offset); + } + } + { + pm["Name"] = "faceoffsets"; + Tag t("DataArray", pm, os); + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << cell_foffsets[c] << ' '; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + { + pm["type"] = "UInt8"; + pm["Name"] = "types"; + Tag t("DataArray", pm, os); + const int num_per_line = 10; + for (int c = 0; c < num_cells; ++c) { + if (c % num_per_line == 0) { + Tag::indent(os); + } + os << "42 "; + if (c % num_per_line == num_per_line - 1 + || c == num_cells - 1) { + os << '\n'; + } + } + } + } + { + pm.clear(); + if (data.find("saturation") != data.end()) { + pm["Scalars"] = "saturation"; + } else if (data.find("pressure") != data.end()) { + pm["Scalars"] = "pressure"; + } + Tag celldatatag("CellData", pm, os); + pm.clear(); + pm["NumberOfComponents"] = "1"; + pm["format"] = "ascii"; + pm["type"] = "Float64"; + for (DataMap::const_iterator dit = data.begin(); dit != data.end(); ++dit) { + pm["Name"] = dit->first; + const std::vector& field = *(dit->second); + const int num_comps = field.size()/grid.number_of_cells; + pm["NumberOfComponents"] = boost::lexical_cast(num_comps); + Tag ptag("DataArray", pm, os); + const int num_per_line = num_comps == 1 ? 5 : num_comps; + for (int item = 0; item < num_cells*num_comps; ++item) { + if (item % num_per_line == 0) { + Tag::indent(os); + } double value = field[item]; if (std::fabs(value) < std::numeric_limits::min()) { // Avoiding denormal numbers to work around // bug in Paraview. value = 0.0; } - os << value << ' '; - if (item % num_per_line == num_per_line - 1 - || item == num_cells - 1) { - os << '\n'; - } - } - } - } + os << value << ' '; + if (item % num_per_line == num_per_line - 1 + || item == num_cells - 1) { + os << '\n'; + } + } + } + } } } // namespace Opm diff --git a/opm/core/utility/writeVtkData.hpp b/opm/core/utility/writeVtkData.hpp index 200b2b26..96cdfac2 100644 --- a/opm/core/utility/writeVtkData.hpp +++ b/opm/core/utility/writeVtkData.hpp @@ -26,25 +26,23 @@ #include #include #include +#include struct UnstructuredGrid; namespace Opm { - /// Intended to map strings (giving the output field names) to data. - typedef std::map*> DataMap; - /// Vtk output for cartesian grids. void writeVtkData(const std::tr1::array& dims, - const std::tr1::array& cell_size, - const DataMap& data, - std::ostream& os); + const std::tr1::array& cell_size, + const DataMap& data, + std::ostream& os); /// Vtk output for general grids. void writeVtkData(const UnstructuredGrid& grid, - const DataMap& data, - std::ostream& os); + const DataMap& data, + std::ostream& os); } // namespace Opm #endif // OPM_WRITEVTKDATA_HEADER_INCLUDED diff --git a/tests/Makefile.am b/tests/Makefile.am index 6e69ae5c..bcf9decb 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -1,8 +1,9 @@ AM_CPPFLAGS = \ -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 @@ -77,3 +78,9 @@ if BUILD_AGMG noinst_PROGRAMS += test_agmg test_agmg_SOURCES = test_agmg.cpp endif + + +if HAVE_ERT +noinst_PROGRAMS += test_ert +test_ert_SOURCES = test_ert.cpp +endif diff --git a/tests/test_ert.cpp b/tests/test_ert.cpp new file mode 100644 index 00000000..282549d1 --- /dev/null +++ b/tests/test_ert.cpp @@ -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 . +*/ + + +#include "config.h" +#include +#include +#include +#include +#include +#include +#include +#include + + +#include + +using namespace std; + +#if 0 +static +void cell_nodes(const UnstructuredGrid * c_grid , int cell , std::vector& 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 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 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 grid; + boost::scoped_ptr 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())); +}