diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake index caa6ca02..75671e11 100644 --- a/CMakeLists_files.cmake +++ b/CMakeLists_files.cmake @@ -165,10 +165,13 @@ list (APPEND TEST_SOURCE_FILES tests/test_param.cpp tests/test_blackoilfluid.cpp tests/test_shadow.cpp - tests/test_units.cpp - tests/test_blackoilstate.cpp - tests/test_wellsmanager.cpp - tests/test_wellcontrols.cpp + tests/test_units.cpp + tests/test_blackoilstate.cpp + tests/test_parser.cpp + tests/test_wellsmanager.cpp + tests/test_wellcontrols.cpp + tests/test_wellsgroup.cpp + tests/test_wellcollection.cpp ) # originally generated with the command: @@ -177,10 +180,13 @@ list (APPEND TEST_DATA_FILES tests/extratestdata.xml tests/testdata.xml tests/liveoil.DATA - tests/wetgas.DATA - tests/testBlackoilState1.DATA - tests/testBlackoilState2.DATA - tests/wells_manager_data.data + tests/wetgas.DATA + tests/testBlackoilState1.DATA + tests/testBlackoilState2.DATA + tests/wells_manager_data.data + tests/wells_manager_data_expanded.data + tests/wells_manager_data_wellSTOP.data + tests/wells_group.data ) # originally generated with the command: diff --git a/cmake/Modules/Findcjson.cmake b/cmake/Modules/Findcjson.cmake index fde5a782..57055739 100644 --- a/cmake/Modules/Findcjson.cmake +++ b/cmake/Modules/Findcjson.cmake @@ -31,13 +31,14 @@ if (CMAKE_SIZEOF_VOID_P) math (EXPR _BITS "8 * ${CMAKE_SIZEOF_VOID_P}") endif (CMAKE_SIZEOF_VOID_P) +string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/?(.*)" "\\1" BUILD_DIR_SUFFIX "${PROJECT_BINARY_DIR}") + find_library (CJSON_LIBRARY NAMES "cjson" HINTS "${CJSON_ROOT}" PATHS "${PROJECT_BINARY_DIR}/../opm-parser" - "${PROJECT_BINARY_DIR}/../opm-parser-build" - "${PROJECT_BINARY_DIR}/../../opm-parser/build" - "${PROJECT_BINARY_DIR}/../../opm-parser/cmake-build" + "${PROJECT_BINARY_DIR}/../opm-parser${BUILD_DIR_SUFFIX}" + "${PROJECT_BINARY_DIR}/../../opm-parser/${BUILD_DIR_SUFFIX}" PATH_SUFFIXES "lib" "lib${_BITS}" "lib/${CMAKE_LIBRARY_ARCHITECTURE}" "opm/json" DOC "Path to cjson library archive/shared object files" diff --git a/cmake/Modules/Findopm-parser.cmake b/cmake/Modules/Findopm-parser.cmake index 3ae6cc0b..2698209f 100644 --- a/cmake/Modules/Findopm-parser.cmake +++ b/cmake/Modules/Findopm-parser.cmake @@ -5,9 +5,9 @@ # # If found, it sets these variables: # -# HAVE_OPM_PARSER Defined if a test program compiled -# OPM_PARSER_INCLUDE_DIRS Header file directories -# OPM_PARSER_LIBRARIES Archives and shared objects +# HAVE_OPM_PARSER Defined if a test program compiled +# OPM_PARSER_INCLUDE_DIRS Header file directories +# OPM_PARSER_LIBRARIES Archives and shared objects include (FindPackageHandleStandardArgs) @@ -31,6 +31,9 @@ if ((NOT OPM_PARSER_ROOT) AND OPM_ROOT) set (OPM_PARSER_ROOT "${OPM_ROOT}/opm-parser") endif () +# Detect the build dir suffix or subdirectory +string(REGEX REPLACE "${PROJECT_SOURCE_DIR}/?(.*)" "\\1" BUILD_DIR_SUFFIX "${PROJECT_BINARY_DIR}") + # if a root is specified, then don't search in system directories # or in relative directories to this one if (OPM_PARSER_ROOT) @@ -43,9 +46,8 @@ else () "${PROJECT_SOURCE_DIR}/../opm-parser") set (_opm_parser_build "${PROJECT_BINARY_DIR}/../opm-parser" - "${PROJECT_BINARY_DIR}/../opm-parser-build" - "${PROJECT_BINARY_DIR}/../../opm-parser/build" - "${PROJECT_BINARY_DIR}/../../opm-parser/cmake-build") + "${PROJECT_BINARY_DIR}/../opm-parser${BUILD_DIR_SUFFIX}" + "${PROJECT_BINARY_DIR}/../../opm-parser/${BUILD_DIR_SUFFIX}") endif () # use this header as signature @@ -61,18 +63,18 @@ find_path (OPM_PARSER_INCLUDE_DIR # then it is probably a build directory; read the CMake cache of # opm-parser to figure out where the source directory is if ((NOT OPM_PARSER_INCLUDE_DIR) AND - (OPM_PARSER_ROOT AND (EXISTS "${OPM_PARSER_ROOT}/CMakeCache.txt"))) + (OPM_PARSER_ROOT AND (EXISTS "${OPM_PARSER_ROOT}/CMakeCache.txt"))) set (_regex "^OPMParser_SOURCE_DIR:STATIC=\(.*\)$") file (STRINGS - "${OPM_PARSER_ROOT}/CMakeCache.txt" - _cache_entry - REGEX "${_regex}") + "${OPM_PARSER_ROOT}/CMakeCache.txt" + _cache_entry + REGEX "${_regex}") string(REGEX REPLACE "${_regex}" "\\1" - OPM_PARSER_INCLUDE_DIR - "${_cache_entry}") + OPM_PARSER_INCLUDE_DIR + "${_cache_entry}") if (OPM_PARSER_INCLUDE_DIR) - set (OPM_PARSER_INCLUDE_DIR "${OPM_PARSER_INCLUDE_DIR}" - CACHE PATH "Path to OPM parser header files" FORCE) + set (OPM_PARSER_INCLUDE_DIR "${OPM_PARSER_INCLUDE_DIR}" + CACHE PATH "Path to OPM parser header files" FORCE) endif () endif () @@ -109,14 +111,16 @@ endif () # get the prerequisite Boost libraries if (NOT Boost_FOUND) find_package(Boost 1.44.0 - COMPONENTS filesystem date_time system unit_test_framework REQUIRED ${OPM_PARSER_QUIET}) + COMPONENTS filesystem date_time system unit_test_framework REQUIRED ${OPM_PARSER_QUIET}) endif () # setup list of all required libraries to link with opm-parser. notice that # we use the plural form to get *all* the libraries needed by cjson set (OPM_PARSER_INCLUDE_DIRS ${OPM_PARSER_INCLUDE_DIR} + ${CJSON_INCLUDE_DIRS} ${Boost_INCLUDE_DIRS}) + set (OPM_PARSER_LIBRARIES ${OPM_PARSER_LIBRARY} ${OPM_JSON_LIBRARY} @@ -126,7 +130,7 @@ set (OPM_PARSER_LIBRARIES # see if we can compile a minimum example # CMake logical test doesn't handle lists (sic) if (NOT (OPM_PARSER_INCLUDE_DIR MATCHES "-NOTFOUND" - OR OPM_PARSER_LIBRARIES MATCHES "-NOTFOUND")) + OR OPM_PARSER_LIBRARIES MATCHES "-NOTFOUND")) include (CMakePushCheckState) include (CheckCSourceCompiles) cmake_push_check_state () diff --git a/cmake/Modules/OpmFind.cmake b/cmake/Modules/OpmFind.cmake index 58501ba9..592a9830 100644 --- a/cmake/Modules/OpmFind.cmake +++ b/cmake/Modules/OpmFind.cmake @@ -59,6 +59,7 @@ set (_opm_proj_exemptions dune-istl dune-grid dune-geometry + opm-parser ) # although a DUNE module, it is delivered in the OPM suite diff --git a/cmake/Modules/opm-core-prereqs.cmake b/cmake/Modules/opm-core-prereqs.cmake index 0c68ec60..41419b15 100644 --- a/cmake/Modules/opm-core-prereqs.cmake +++ b/cmake/Modules/opm-core-prereqs.cmake @@ -3,31 +3,33 @@ # defines that must be present in config.h for our headers set (opm-core_CONFIG_VAR - HAVE_ERT - HAVE_SUITESPARSE_UMFPACK_H - ) + HAVE_ERT + HAVE_SUITESPARSE_UMFPACK_H + ) # dependencies set (opm-core_DEPS - # compile with C99 support if available - "C99" - # compile with C++0x/11 support if available - "CXX11Features REQUIRED" - # various runtime library enhancements - "Boost 1.39.0 - COMPONENTS date_time filesystem system unit_test_framework REQUIRED" - # matrix library - "BLAS REQUIRED" - "LAPACK REQUIRED" - # Tim Davis' SuiteSparse archive - "SuiteSparse COMPONENTS umfpack" - # solver - "SuperLU" - # xml processing (for config parsing) - "TinyXML" - # Ensembles-based Reservoir Tools (ERT) - "ERT" - # DUNE dependency - "dune-common" - "dune-istl" - ) + # compile with C99 support if available + "C99" + # compile with C++0x/11 support if available + "CXX11Features REQUIRED" + # various runtime library enhancements + "Boost 1.44.0 + COMPONENTS date_time filesystem system unit_test_framework REQUIRED" + # matrix library + "BLAS REQUIRED" + "LAPACK REQUIRED" + # Tim Davis' SuiteSparse archive + "SuiteSparse COMPONENTS umfpack" + # solver + "SuperLU" + # xml processing (for config parsing) + "TinyXML" + #Parser library + "opm-parser REQUIRED" + # Ensembles-based Reservoir Tools (ERT) + "ERT" + # DUNE dependency + "dune-common" + "dune-istl" + ) diff --git a/cmake/Modules/opm-parser-prereqs.cmake b/cmake/Modules/opm-parser-prereqs.cmake index 80c95136..521568b7 100644 --- a/cmake/Modules/opm-parser-prereqs.cmake +++ b/cmake/Modules/opm-parser-prereqs.cmake @@ -13,6 +13,6 @@ set (opm-parser_DEPS # compile with C++0x/11 support if available "CXX11Features REQUIRED" # various runtime library enhancements - "Boost 1.39.0 - COMPONENTS date_time filesystem system unit_test_framework REQUIRED" + "Boost 1.44.0 COMPONENTS date_time filesystem system unit_test_framework REQUIRED" + "cJSON" ) diff --git a/dune.module b/dune.module index 02dc99c3..b3b7ad9c 100644 --- a/dune.module +++ b/dune.module @@ -4,4 +4,4 @@ Description: Open Porous Media Initiative Core Library Version: 1.1 Label: 2013.10 Maintainer: atgeirr@sintef.no -Depends: dune-common (>= 2.2) dune-istl (>= 2.2) +Depends: dune-common (>= 2.2) dune-istl (>= 2.2) opm-parser diff --git a/opm/core/grid/GridManager.cpp b/opm/core/grid/GridManager.cpp index f557f0ca..392de07c 100644 --- a/opm/core/grid/GridManager.cpp +++ b/opm/core/grid/GridManager.cpp @@ -23,16 +23,13 @@ #include #include #include + +#include #include #include - - namespace Opm { - - - /// Construct a 3d corner-point grid from a deck. GridManager::GridManager(const Opm::EclipseGridParser& deck) { @@ -58,6 +55,30 @@ namespace Opm } } + /// Construct a 3d corner-point grid from a deck. + GridManager::GridManager(Opm::DeckConstPtr newParserDeck) + { + // We accept two different ways to specify the grid. + // 1. Corner point format. + // 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. + // Requires DXV, DYV, DZV, optionally DEPTHZ or TOPS. + // For this format, we will verify that ZCORN, COORDS + // and ACTNUM are not present. + // Note that for TOPS, we only allow a uniform vector of values. + + if (newParserDeck->hasKeyword("ZCORN") && newParserDeck->hasKeyword("COORD")) { + initFromDeckCornerpoint(newParserDeck); + } else if (newParserDeck->hasKeyword("DXV") && newParserDeck->hasKeyword("DYV") && newParserDeck->hasKeyword("DZV")) { + initFromDeckTensorgrid(newParserDeck); + } else { + OPM_THROW(std::runtime_error, "Could not initialize grid from deck. " + "Need either ZCORN + COORD or DXV + DYV + DZV keywords."); + } + } @@ -152,6 +173,71 @@ namespace Opm } } + // Construct corner-point grid from deck. + void GridManager::initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck) + { + // Extract data from deck. + // Collect in input struct for preprocessing. + struct grdecl grdecl; + createGrdecl(newParserDeck, grdecl); + + // Process grid. + ug_ = create_grid_cornerpoint(&grdecl, 0.0); + if (!ug_) { + OPM_THROW(std::runtime_error, "Failed to construct grid."); + } + } + + void GridManager::createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl) + { + // Extract data from deck. + const std::vector& zcorn = newParserDeck->getKeyword("ZCORN")->getSIDoubleData(); + const std::vector& coord = newParserDeck->getKeyword("COORD")->getSIDoubleData(); + const int* actnum = NULL; + if (newParserDeck->hasKeyword("ACTNUM")) { + actnum = &(newParserDeck->getKeyword("ACTNUM")->getIntData()[0]); + } + + std::array dims; + if (newParserDeck->hasKeyword("DIMENS")) { + Opm::DeckKeywordConstPtr dimensKeyword = newParserDeck->getKeyword("DIMENS"); + dims[0] = dimensKeyword->getRecord(0)->getItem(0)->getInt(0); + dims[1] = dimensKeyword->getRecord(0)->getItem(1)->getInt(0); + dims[2] = dimensKeyword->getRecord(0)->getItem(2)->getInt(0); + } else if (newParserDeck->hasKeyword("SPECGRID")) { + Opm::DeckKeywordConstPtr specgridKeyword = newParserDeck->getKeyword("SPECGRID"); + dims[0] = specgridKeyword->getRecord(0)->getItem(0)->getInt(0); + dims[1] = specgridKeyword->getRecord(0)->getItem(1)->getInt(0); + dims[2] = specgridKeyword->getRecord(0)->getItem(2)->getInt(0); + } else { + OPM_THROW(std::runtime_error, "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 (newParserDeck->hasKeyword("MAPAXES")) { + Opm::DeckKeywordConstPtr mapaxesKeyword = newParserDeck->getKeyword("MAPAXES"); + Opm::DeckRecordConstPtr mapaxesRecord = mapaxesKeyword->getRecord(0); + + // memleak alert: here we need to make sure that C code + // can properly take ownership of the grdecl.mapaxes + // object. if it is not freed, it will result in a + // memleak... + double *cWtfMapaxes = static_cast(malloc(sizeof(double)*mapaxesRecord->size())); + for (unsigned i = 0; i < mapaxesRecord->size(); ++i) + cWtfMapaxes[i] = mapaxesRecord->getItem(i)->getSIDouble(0); + grdecl.mapaxes = cWtfMapaxes; + } else + grdecl.mapaxes = NULL; + + } namespace { @@ -230,6 +316,71 @@ namespace Opm } } + void GridManager::initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck) + { + // Extract logical cartesian size. + std::array dims; + if (newParserDeck->hasKeyword("DIMENS")) { + Opm::DeckKeywordConstPtr dimensKeyword = newParserDeck->getKeyword("DIMENS"); + dims[0] = dimensKeyword->getRecord(0)->getItem(0)->getInt(0); + dims[1] = dimensKeyword->getRecord(0)->getItem(1)->getInt(0); + dims[2] = dimensKeyword->getRecord(0)->getItem(2)->getInt(0); + } else if (newParserDeck->hasKeyword("SPECGRID")) { + Opm::DeckKeywordConstPtr specgridKeyword = newParserDeck->getKeyword("SPECGRID"); + dims[0] = specgridKeyword->getRecord(0)->getItem(0)->getInt(0); + dims[1] = specgridKeyword->getRecord(0)->getItem(1)->getInt(0); + dims[2] = specgridKeyword->getRecord(0)->getItem(2)->getInt(0); + } else { + OPM_THROW(std::runtime_error, "Deck must have either DIMENS or SPECGRID."); + } + + // Extract coordinates (or offsets from top, in case of z). + const std::vector& dxv = newParserDeck->getKeyword("DXV")->getSIDoubleData(); + const std::vector& dyv = newParserDeck->getKeyword("DYV")->getSIDoubleData(); + const std::vector& dzv = newParserDeck->getKeyword("DZV")->getSIDoubleData(); + std::vector x = coordsFromDeltas(dxv); + std::vector y = coordsFromDeltas(dyv); + std::vector z = coordsFromDeltas(dzv); + + // Check that number of cells given are consistent with DIMENS/SPECGRID. + if (dims[0] != int(dxv.size())) { + OPM_THROW(std::runtime_error, "Number of DXV data points do not match DIMENS or SPECGRID."); + } + if (dims[1] != int(dyv.size())) { + OPM_THROW(std::runtime_error, "Number of DYV data points do not match DIMENS or SPECGRID."); + } + if (dims[2] != int(dzv.size())) { + OPM_THROW(std::runtime_error, "Number of DZV data points do not match DIMENS or SPECGRID."); + } + + // Extract top corner depths, if available. + const double* top_depths = 0; + std::vector top_depths_vec; + if (newParserDeck->hasKeyword("DEPTHZ")) { + const std::vector& depthz = newParserDeck->getKeyword("DEPTHZ")->getSIDoubleData(); + if (depthz.size() != x.size()*y.size()) { + OPM_THROW(std::runtime_error, "Incorrect size of DEPTHZ: " << depthz.size()); + } + top_depths = &depthz[0]; + } else if (newParserDeck->hasKeyword("TOPS")) { + // We only support constant values for TOPS. + // It is not 100% clear how we best can deal with + // varying TOPS (stair-stepping grid, or not). + const std::vector& tops = newParserDeck->getKeyword("TOPS")->getSIDoubleData(); + if (std::count(tops.begin(), tops.end(), tops[0]) != int(tops.size())) { + OPM_THROW(std::runtime_error, "We do not support nonuniform TOPS, please use ZCORN/COORDS instead."); + } + top_depths_vec.resize(x.size()*y.size(), tops[0]); + top_depths = &top_depths_vec[0]; + } + + // Construct grid. + ug_ = create_grid_tensor3d(dxv.size(), dyv.size(), dzv.size(), + &x[0], &y[0], &z[0], top_depths); + if (!ug_) { + OPM_THROW(std::runtime_error, "Failed to construct grid."); + } + } } // namespace Opm diff --git a/opm/core/grid/GridManager.hpp b/opm/core/grid/GridManager.hpp index 03666377..030d5555 100644 --- a/opm/core/grid/GridManager.hpp +++ b/opm/core/grid/GridManager.hpp @@ -20,10 +20,12 @@ #ifndef OPM_GRIDMANAGER_HEADER_INCLUDED #define OPM_GRIDMANAGER_HEADER_INCLUDED +#include + #include struct UnstructuredGrid; - +struct grdecl; namespace Opm { @@ -44,6 +46,9 @@ namespace Opm /// Construct a 3d corner-point grid or tensor grid from a deck. GridManager(const Opm::EclipseGridParser& deck); + /// Construct a 3d corner-point grid or tensor grid from a deck. + GridManager(Opm::DeckConstPtr newParserDeck); + /// Construct a 2d cartesian grid with cells of unit size. GridManager(int nx, int ny); @@ -72,6 +77,8 @@ namespace Opm /// to make it clear that we are returning a C-compatible struct. const UnstructuredGrid* c_grid() const; + static void createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl); + private: // Disable copying and assignment. GridManager(const GridManager& other); @@ -79,8 +86,10 @@ namespace Opm // Construct corner-point grid from deck. void initFromDeckCornerpoint(const Opm::EclipseGridParser& deck); + void initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck); // Construct tensor grid from deck. void initFromDeckTensorgrid(const Opm::EclipseGridParser& deck); + void initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck); // The managed UnstructuredGrid. UnstructuredGrid* ug_; diff --git a/opm/core/props/BlackoilPropertiesFromDeck.cpp b/opm/core/props/BlackoilPropertiesFromDeck.cpp index caa27f1b..b901ccfb 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.cpp @@ -23,7 +23,6 @@ namespace Opm { - BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, bool init_rock) @@ -43,6 +42,25 @@ namespace Opm } } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + bool init_rock) + { + if (init_rock){ + rock_.init(newParserDeck, grid); + } + pvt_.init(newParserDeck, /*numSamples=*/200); + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, /*numSamples=*/200); + + if (pvt_.numPhases() != satprops_->numPhases()) { + OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); + } + } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck, const UnstructuredGrid& grid, const parameter::ParameterGroup& param, @@ -107,6 +125,70 @@ namespace Opm } } + BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param, + bool init_rock) + { + if(init_rock){ + rock_.init(newParserDeck, grid); + } + + const int pvt_samples = param.getDefault("pvt_tab_size", 200); + pvt_.init(newParserDeck, pvt_samples); + + // Unfortunate lack of pointer smartness here... + const int sat_samples = param.getDefault("sat_tab_size", 200); + std::string threephase_model = param.getDefault("threephase_model", "simple"); + if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "simple") { + OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only."); + } + if (sat_samples > 1) { + if (threephase_model == "stone2") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "simple") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else { + OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); + } + } else { + if (threephase_model == "stone2") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "simple") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(newParserDeck, grid, sat_samples); + } else { + OPM_THROW(std::runtime_error, "Unknown threephase_model: " << threephase_model); + } + } + + if (pvt_.numPhases() != satprops_->numPhases()) { + OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data (" + << pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ")."); + } + } + BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck() { } diff --git a/opm/core/props/BlackoilPropertiesFromDeck.hpp b/opm/core/props/BlackoilPropertiesFromDeck.hpp index 49072b29..5c74a76a 100644 --- a/opm/core/props/BlackoilPropertiesFromDeck.hpp +++ b/opm/core/props/BlackoilPropertiesFromDeck.hpp @@ -27,6 +27,9 @@ #include #include #include + +#include + #include struct UnstructuredGrid; @@ -44,7 +47,15 @@ namespace Opm /// \param[in] grid Grid to which property object applies, needed for the /// mapping from cell indices (typically from a processed grid) /// to logical cartesian indices consistent with the deck. - BlackoilPropertiesFromDeck(const EclipseGridParser& deck, + BlackoilPropertiesFromDeck(const EclipseGridParser &deck, + const UnstructuredGrid& grid, bool init_rock=true ); + + /// Initialize from deck and grid. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, const UnstructuredGrid& grid, bool init_rock=true ); /// Initialize from deck, grid and parameters. @@ -63,6 +74,22 @@ namespace Opm const parameter::ParameterGroup& param, bool init_rock=true); + /// Initialize from deck, grid and parameters. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + /// \param[in] param Parameters. Accepted parameters include: + /// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables. + /// sat_tab_size (200) number of uniform sample points for saturation tables. + /// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2"). + /// For both size parameters, a 0 or negative value indicates that no spline fitting is to + /// be done, and the input fluid data used directly for linear interpolation. + BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const parameter::ParameterGroup& param, + bool init_rock=true); + /// Destructor. virtual ~BlackoilPropertiesFromDeck(); diff --git a/opm/core/props/phaseUsageFromDeck.hpp b/opm/core/props/phaseUsageFromDeck.hpp index 01f28ee1..6ee66075 100644 --- a/opm/core/props/phaseUsageFromDeck.hpp +++ b/opm/core/props/phaseUsageFromDeck.hpp @@ -24,10 +24,49 @@ #include #include +#include +#include + namespace Opm { + /// Looks at presence of WATER, OIL and GAS keywords in state object + /// to determine active phases. + inline PhaseUsage phaseUsageFromDeck(Opm::EclipseStateConstPtr eclipseState) + { + PhaseUsage pu; + std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0); + + // Discover phase usage. + if (eclipseState->hasPhase(Phase::PhaseEnum::WATER)) { + pu.phase_used[BlackoilPhases::Aqua] = 1; + } + if (eclipseState->hasPhase(Phase::PhaseEnum::OIL)) { + pu.phase_used[BlackoilPhases::Liquid] = 1; + } + if (eclipseState->hasPhase(Phase::PhaseEnum::GAS)) { + pu.phase_used[BlackoilPhases::Vapour] = 1; + } + pu.num_phases = 0; + for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) { + pu.phase_pos[i] = pu.num_phases; + pu.num_phases += pu.phase_used[i]; + } + + // Only 2 or 3 phase systems handled. + if (pu.num_phases < 2 || pu.num_phases > 3) { + OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases."); + } + + // We need oil systems, since we do not support the keywords needed for + // water-gas systems. + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems."); + } + + return pu; + } /// Looks at presence of WATER, OIL and GAS keywords in deck /// to determine active phases. @@ -66,6 +105,43 @@ namespace Opm return pu; } + /// Looks at presence of WATER, OIL and GAS keywords in deck + /// to determine active phases. + inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr newParserDeck) + { + PhaseUsage pu; + std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0); + + // Discover phase usage. + if (newParserDeck->hasKeyword("WATER")) { + pu.phase_used[BlackoilPhases::Aqua] = 1; + } + if (newParserDeck->hasKeyword("OIL")) { + pu.phase_used[BlackoilPhases::Liquid] = 1; + } + if (newParserDeck->hasKeyword("GAS")) { + pu.phase_used[BlackoilPhases::Vapour] = 1; + } + pu.num_phases = 0; + for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) { + pu.phase_pos[i] = pu.num_phases; + pu.num_phases += pu.phase_used[i]; + } + + // Only 2 or 3 phase systems handled. + if (pu.num_phases < 2 || pu.num_phases > 3) { + OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases."); + } + + // We need oil systems, since we do not support the keywords needed for + // water-gas systems. + if (!pu.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems."); + } + + return pu; + } + } #endif // OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED diff --git a/opm/core/props/pvt/BlackoilPvtProperties.cpp b/opm/core/props/pvt/BlackoilPvtProperties.cpp index 217107cb..0fd7dda4 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.cpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.cpp @@ -32,6 +32,9 @@ #include #include +#include +#include +#include namespace Opm { @@ -40,7 +43,6 @@ namespace Opm { } - void BlackoilPvtProperties::init(const EclipseGridParser& deck, const int samples) { // If we need multiple regions, this class and the SinglePvt* classes must change. @@ -113,6 +115,83 @@ namespace Opm } } + void BlackoilPvtProperties::init(Opm::DeckConstPtr newParserDeck, int samples) + { + // If we need multiple regions, this class and the SinglePvt* classes must change. + region_number_ = 0; + + phase_usage_ = phaseUsageFromDeck(newParserDeck); + + // Surface densities. Accounting for different orders in eclipse and our code. + if (newParserDeck->hasKeyword("DENSITY")) { + Opm::DeckKeywordConstPtr densityKeyword = newParserDeck->getKeyword("DENSITY"); + const std::vector& d = densityKeyword->getRecord(region_number_)->getItem(0)->getSIDoubleData(); + enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 }; + if (phase_usage_.phase_used[Aqua]) { + densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water]; + } + if (phase_usage_.phase_used[Vapour]) { + densities_[phase_usage_.phase_pos[Vapour]] = d[ECL_gas]; + } + if (phase_usage_.phase_used[Liquid]) { + densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil]; + } + } else { + OPM_THROW(std::runtime_error, "Input is missing DENSITY\n"); + } + + // Set the properties. + props_.resize(phase_usage_.num_phases); + // Water PVT + if (phase_usage_.phase_used[Aqua]) { + if (newParserDeck->hasKeyword("PVTW")) { + Opm::PvtwTable pvtwTable(newParserDeck->getKeyword("PVTW")); + + props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable)); + } else { + // Eclipse 100 default. + props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise)); + } + } + // Oil PVT + if (phase_usage_.phase_used[Liquid]) { + if (newParserDeck->hasKeyword("PVDO")) { + Opm::PvdoTable pvdoTable(newParserDeck->getKeyword("PVDO"), region_number_); + + props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples)); + } else if (newParserDeck->hasKeyword("PVTO")) { + Opm::PvtoTable pvtoTable(newParserDeck->getKeyword("PVTO"), /*tableIdx=*/0); + + props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable)); + } else if (newParserDeck->hasKeyword("PVCDO")) { + Opm::PvdcoTable pvcdoTable(newParserDeck->getKeyword("PVCDO")); + + props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable)); + } else { + OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n"); + } + } + // Gas PVT + if (phase_usage_.phase_used[Vapour]) { + if (newParserDeck->hasKeyword("PVDG")) { + Opm::PvdgTable pvdgTable(newParserDeck->getKeyword("PVDG"), region_number_); + + props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples)); + } else if (newParserDeck->hasKeyword("PVTG")) { + Opm::PvtgTable pvtgTable(newParserDeck->getKeyword("PVTG"), /*tableIdx=*/0); + + props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable)); + } else { + OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n"); + } + } + + // Must inform pvt property objects of phase structure. + for (int i = 0; i < phase_usage_.num_phases; ++i) { + props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos); + } + } + const double* BlackoilPvtProperties::surfaceDensities() const { return densities_; diff --git a/opm/core/props/pvt/BlackoilPvtProperties.hpp b/opm/core/props/pvt/BlackoilPvtProperties.hpp index 4ab60921..9b466b2d 100644 --- a/opm/core/props/pvt/BlackoilPvtProperties.hpp +++ b/opm/core/props/pvt/BlackoilPvtProperties.hpp @@ -20,11 +20,12 @@ #ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED #define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED - - #include #include #include + +#include + #include #include @@ -55,6 +56,11 @@ namespace Opm /// data without fitting a spline. void init(const EclipseGridParser& deck, const int samples); + + /// Initialize from deck. + /// \param deck An input deck from the opm-parser module. + void init(Opm::DeckConstPtr deck, int samples); + /// \return Object describing the active phases. PhaseUsage phaseUsage() const; diff --git a/opm/core/props/pvt/SinglePvtConstCompr.hpp b/opm/core/props/pvt/SinglePvtConstCompr.hpp index e97a7242..d6c4d088 100644 --- a/opm/core/props/pvt/SinglePvtConstCompr.hpp +++ b/opm/core/props/pvt/SinglePvtConstCompr.hpp @@ -20,9 +20,12 @@ #ifndef OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED #define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED - #include #include + +#include +#include + #include #include @@ -55,6 +58,32 @@ namespace Opm visc_comp_ = pvtw[region_number][4]; } + SinglePvtConstCompr(const Opm::PvtwTable &pvtwTable) + { + if (pvtwTable.numRows() != 1) + OPM_THROW(std::runtime_error, + "The table specified by the PVTW keyword is required" + "to exhibit exactly one row."); + ref_press_ = pvtwTable.getPressureColumn()[0]; + ref_B_ = pvtwTable.getFormationFactorColumn()[0]; + comp_ = pvtwTable.getCompressibilityColumn()[0]; + viscosity_ = pvtwTable.getViscosityColumn()[0]; + visc_comp_ = pvtwTable.getViscosibilityColumn()[0]; + } + + SinglePvtConstCompr(const Opm::PvdcoTable &pvdcoTable) + { + if (pvdcoTable.numRows() != 1) + OPM_THROW(std::runtime_error, + "The table specified by the PVDCO keyword is required" + "to exhibit exactly one row."); + ref_press_ = pvdcoTable.getPressureColumn()[0]; + ref_B_ = pvdcoTable.getFormationFactorColumn()[0]; + comp_ = pvdcoTable.getCompressibilityColumn()[0]; + viscosity_ = pvdcoTable.getViscosityColumn()[0]; + visc_comp_ = pvdcoTable.getViscosibilityColumn()[0]; + } + SinglePvtConstCompr(double visc) : ref_press_(0.0), ref_B_(1.0), diff --git a/opm/core/props/pvt/SinglePvtDead.cpp b/opm/core/props/pvt/SinglePvtDead.cpp index 2db23dbe..fb11d2de 100644 --- a/opm/core/props/pvt/SinglePvtDead.cpp +++ b/opm/core/props/pvt/SinglePvtDead.cpp @@ -62,6 +62,30 @@ namespace Opm // << "\n\nvisc\n\n" << viscosity_ << std::endl; } + /// Constructor + SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable) + { + // Copy data + const std::vector& press = pvdoTable.getPressureColumn(); + const std::vector& b = pvdoTable.getFormationFactorColumn(); + const std::vector& visc = pvdoTable.getViscosityColumn(); + + const int sz = b.size(); + std::vector bInv(sz); + for (int i = 0; i < sz; ++i) { + bInv[i] = 1.0 / b[i]; + } + b_ = NonuniformTableLinear(press, bInv); + viscosity_ = NonuniformTableLinear(press, visc); + + // Dumping the created tables. +// static int count = 0; +// std::ofstream os((std::string("dump-") + boost::lexical_cast(count++)).c_str()); +// os.precision(15); +// os << "1/B\n\n" << one_over_B_ +// << "\n\nvisc\n\n" << viscosity_ << std::endl; + } + // Destructor SinglePvtDead::~SinglePvtDead() { diff --git a/opm/core/props/pvt/SinglePvtDead.hpp b/opm/core/props/pvt/SinglePvtDead.hpp index 2c40a951..703da0b8 100644 --- a/opm/core/props/pvt/SinglePvtDead.hpp +++ b/opm/core/props/pvt/SinglePvtDead.hpp @@ -23,6 +23,9 @@ #include #include + +#include + #include namespace Opm @@ -40,6 +43,7 @@ namespace Opm public: typedef std::vector > > table_t; SinglePvtDead(const table_t& pvd_table); + SinglePvtDead(const Opm::PvdoTable &pvdoTable); virtual ~SinglePvtDead(); /// Viscosity as a function of p and z. diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.cpp b/opm/core/props/pvt/SinglePvtDeadSpline.cpp index a0873e0b..0dc3b560 100644 --- a/opm/core/props/pvt/SinglePvtDeadSpline.cpp +++ b/opm/core/props/pvt/SinglePvtDeadSpline.cpp @@ -20,6 +20,9 @@ #include "config.h" #include #include + +#include + #include // Extra includes for debug dumping of tables. @@ -29,7 +32,6 @@ namespace Opm { - //------------------------------------------------------------------------ // Member functions //------------------------------------------------------------------------- @@ -54,13 +56,42 @@ namespace Opm } buildUniformMonotoneTable(press, B_inv, samples, b_); buildUniformMonotoneTable(press, visc, samples, viscosity_); + } - // Dumping the created tables. -// static int count = 0; -// std::ofstream os((std::string("dump-") + boost::lexical_cast(count++)).c_str()); -// os.precision(15); -// os << "1/B\n\n" << one_over_B_ -// << "\n\nvisc\n\n" << viscosity_ << std::endl; + SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples) + { + int numRows = pvdoTable.numRows(); + + // Copy data + const std::vector &press = pvdoTable.getPressureColumn(); + const std::vector &B = pvdoTable.getFormationFactorColumn(); + const std::vector &visc = pvdoTable.getViscosityColumn(); + + std::vector B_inv(numRows); + for (int i = 0; i < numRows; ++i) { + B_inv[i] = 1.0 / B[i]; + } + + buildUniformMonotoneTable(press, B_inv, samples, b_); + buildUniformMonotoneTable(press, visc, samples, viscosity_); + } + + SinglePvtDeadSpline::SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples) + { + int numRows = pvdgTable.numRows(); + + // Copy data + const std::vector &press = pvdgTable.getPressureColumn(); + const std::vector &B = pvdgTable.getFormationFactorColumn(); + const std::vector &visc = pvdgTable.getViscosityColumn(); + + std::vector B_inv(numRows); + for (int i = 0; i < numRows; ++i) { + B_inv[i] = 1.0 / B[i]; + } + + buildUniformMonotoneTable(press, B_inv, samples, b_); + buildUniformMonotoneTable(press, visc, samples, viscosity_); } // Destructor diff --git a/opm/core/props/pvt/SinglePvtDeadSpline.hpp b/opm/core/props/pvt/SinglePvtDeadSpline.hpp index 337c8a80..65bd2561 100644 --- a/opm/core/props/pvt/SinglePvtDeadSpline.hpp +++ b/opm/core/props/pvt/SinglePvtDeadSpline.hpp @@ -20,9 +20,12 @@ #ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED #define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED - #include #include + +#include +#include + #include namespace Opm @@ -41,6 +44,8 @@ namespace Opm typedef std::vector > > table_t; SinglePvtDeadSpline(const table_t& pvd_table, const int samples); + SinglePvtDeadSpline(const Opm::PvdoTable &pvdoTable, int samples); + SinglePvtDeadSpline(const Opm::PvdgTable &pvdgTable, int samples); virtual ~SinglePvtDeadSpline(); /// Viscosity as a function of p and z. diff --git a/opm/core/props/pvt/SinglePvtLiveGas.cpp b/opm/core/props/pvt/SinglePvtLiveGas.cpp index 6a43f475..43debe7b 100644 --- a/opm/core/props/pvt/SinglePvtLiveGas.cpp +++ b/opm/core/props/pvt/SinglePvtLiveGas.cpp @@ -32,8 +32,8 @@ #include #include #include -#include +#include namespace Opm { @@ -41,6 +41,7 @@ namespace Opm using Opm::linearInterpolation; using Opm::linearInterpolationDerivative; + //------------------------------------------------------------------------ // Member functions //------------------------------------------------------------------------- @@ -81,6 +82,27 @@ namespace Opm } } + SinglePvtLiveGas::SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable) + { + // GAS, PVTG + saturated_gas_table_.resize(4); + saturated_gas_table_[0] = pvtgTable.getOuterTable()->getPressureColumn(); + saturated_gas_table_[1] = pvtgTable.getOuterTable()->getGasFormationFactorColumn(); + saturated_gas_table_[2] = pvtgTable.getOuterTable()->getGasViscosityColumn(); + saturated_gas_table_[3] = pvtgTable.getOuterTable()->getOilSolubilityColumn(); + + int sz = pvtgTable.getOuterTable()->numRows(); + undersat_gas_tables_.resize(sz); + for (int i=0; igetGasViscosityColumn(); + } + } + // Destructor SinglePvtLiveGas::~SinglePvtLiveGas() { diff --git a/opm/core/props/pvt/SinglePvtLiveGas.hpp b/opm/core/props/pvt/SinglePvtLiveGas.hpp index de595393..c09f8ffc 100644 --- a/opm/core/props/pvt/SinglePvtLiveGas.hpp +++ b/opm/core/props/pvt/SinglePvtLiveGas.hpp @@ -21,6 +21,9 @@ #define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED #include + +#include + #include namespace Opm @@ -38,6 +41,7 @@ namespace Opm typedef std::vector > > table_t; SinglePvtLiveGas(const table_t& pvto); + SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable); virtual ~SinglePvtLiveGas(); /// Viscosity as a function of p and z. diff --git a/opm/core/props/pvt/SinglePvtLiveOil.cpp b/opm/core/props/pvt/SinglePvtLiveOil.cpp index d18252ee..66f1036e 100644 --- a/opm/core/props/pvt/SinglePvtLiveOil.cpp +++ b/opm/core/props/pvt/SinglePvtLiveOil.cpp @@ -21,8 +21,8 @@ #include #include #include -#include +#include namespace Opm { @@ -31,7 +31,6 @@ namespace Opm using Opm::linearInterpolationDerivative; using Opm::tableIndex; - //------------------------------------------------------------------------ // Member functions //------------------------------------------------------------------------- @@ -102,7 +101,73 @@ namespace Opm undersat_oil_tables_[i][2].push_back(mu); } } + } + SinglePvtLiveOil::SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable) + { + const auto saturatedPvto = pvtoTable.getOuterTable(); + + // OIL, PVTO + saturated_oil_table_.resize(4); + const int sz = saturatedPvto->numRows(); + for (int k=0; k<4; ++k) { + saturated_oil_table_[k].resize(sz); + } + for (int i=0; igetPressureColumn()[i]; // p + saturated_oil_table_[1][i] = 1.0/saturatedPvto->getOilFormationFactorColumn()[i]; // 1/Bo + saturated_oil_table_[2][i] = saturatedPvto->getOilViscosityColumn()[i]; // mu_o + saturated_oil_table_[3][i] = saturatedPvto->getGasSolubilityColumn()[i]; // Rs + } + + undersat_oil_tables_.resize(sz); + for (int i=0; inumRows(); + undersat_oil_tables_[i][0].resize(tsize); + undersat_oil_tables_[i][1].resize(tsize); + undersat_oil_tables_[i][2].resize(tsize); + for (int j=0; jgetPressureColumn()[j]; // p + undersat_oil_tables_[i][1][j] = 1.0/undersaturatedPvto->getOilFormationFactorColumn()[j]; // 1/Bo + undersat_oil_tables_[i][2][j] = undersaturatedPvto->getOilViscosityColumn()[j]; // mu_o + } + } + + // Complete undersaturated tables by extrapolating from existing data + // as is done in Eclipse and Mrst + int iNext = -1; + for (int i=0; i 1) { + continue; + } + // Look ahead for next record containing undersaturated data + if (iNext < i) { + iNext = i+1; + while (iNext > >::size_type sz_t; + for (sz_t j=1; j + +#include + #include namespace Opm @@ -39,6 +41,7 @@ namespace Opm typedef std::vector > > table_t; SinglePvtLiveOil(const table_t& pvto); + SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable); virtual ~SinglePvtLiveOil(); /// Viscosity as a function of p and z. diff --git a/opm/core/props/rock/RockCompressibility.cpp b/opm/core/props/rock/RockCompressibility.cpp index 7b0f69cb..d1729533 100644 --- a/opm/core/props/rock/RockCompressibility.cpp +++ b/opm/core/props/rock/RockCompressibility.cpp @@ -25,6 +25,9 @@ #include #include +#include +#include + #include namespace Opm @@ -65,6 +68,44 @@ namespace Opm } } + RockCompressibility::RockCompressibility(Opm::DeckConstPtr newParserDeck) + : pref_(0.0), + rock_comp_(0.0) + { + if (newParserDeck->hasKeyword("ROCKTAB")) { + Opm::DeckKeywordConstPtr rtKeyword = newParserDeck->getKeyword("ROCKTAB"); + if (rtKeyword->size() != 1) + OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB."); + + // the number of colums of the "ROCKTAB" keyword + // depends on the presence of the "RKTRMDIR" + // keyword. Messy stuff... + bool isDirectional = newParserDeck->hasKeyword("RKTRMDIR"); + if (isDirectional) + { + // well, okay. we don't support non-isotropic + // transmissibility multipliers yet + OPM_THROW(std::runtime_error, "Support for non-isotropic " + "transmissibility multipliers is not implemented yet."); + }; + + Opm::RocktabTable rocktabTable(rtKeyword, isDirectional); + + p_ = rocktabTable.getPressureColumn(); + poromult_ = rocktabTable.getPoreVolumeMultiplierColumn(); + transmult_ = rocktabTable.getTransmissibilityMultiplierColumn(); + } else if (newParserDeck->hasKeyword("ROCK")) { + Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK")); + if (rockTable.numRows() != 1) + OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK."); + + pref_ = rockTable.getPressureColumn()[0]; + rock_comp_ = rockTable.getCompressibilityColumn()[0]; + } else { + std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl; + } + } + bool RockCompressibility::isActive() const { return !p_.empty() || (rock_comp_ != 0.0); diff --git a/opm/core/props/rock/RockCompressibility.hpp b/opm/core/props/rock/RockCompressibility.hpp index 03567cae..bd375e1c 100644 --- a/opm/core/props/rock/RockCompressibility.hpp +++ b/opm/core/props/rock/RockCompressibility.hpp @@ -20,6 +20,8 @@ #ifndef OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED #define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED +#include + #include namespace Opm @@ -35,6 +37,10 @@ namespace Opm /// Looks for the keywords ROCK and ROCKTAB. RockCompressibility(const EclipseGridParser& deck); + /// Construct from input deck. + /// Looks for the keywords ROCK and ROCKTAB. + RockCompressibility(Opm::DeckConstPtr newParserDeck); + /// Construct from parameters. /// Accepts the following parameters (with defaults). /// rock_compressibility_pref (100.0) [given in bar] diff --git a/opm/core/props/rock/RockFromDeck.cpp b/opm/core/props/rock/RockFromDeck.cpp index 73036343..1058fa99 100644 --- a/opm/core/props/rock/RockFromDeck.cpp +++ b/opm/core/props/rock/RockFromDeck.cpp @@ -21,6 +21,9 @@ #include "config.h" #include #include + +#include + #include namespace Opm @@ -37,6 +40,10 @@ namespace Opm PermeabilityKind fillTensor(const EclipseGridParser& parser, std::vector*>& tensor, std::array& kmap); + PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck, + std::vector*>& tensor, + std::array& kmap); + } // anonymous namespace @@ -64,6 +71,15 @@ namespace Opm assignPermeability(deck, grid, perm_threshold); } + void RockFromDeck::init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + assignPorosity(newParserDeck, grid); + permfield_valid_.assign(grid.number_of_cells, false); + const double perm_threshold = 0.0; // Maybe turn into parameter? + assignPermeability(newParserDeck, grid, perm_threshold); + } + void RockFromDeck::assignPorosity(const EclipseGridParser& parser, const UnstructuredGrid& grid) @@ -79,6 +95,21 @@ namespace Opm } } + void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + porosity_.assign(grid.number_of_cells, 1.0); + const int* gc = grid.global_cell; + if (newParserDeck->hasKeyword("PORO")) { + const std::vector& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData(); + for (int c = 0; c < int(porosity_.size()); ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + assert(0 <= c && c < (int) porosity_.size()); + assert(0 <= deck_pos && deck_pos < (int) poro.size()); + porosity_[c] = poro[deck_pos]; + } + } + } void RockFromDeck::assignPermeability(const EclipseGridParser& parser, @@ -135,6 +166,59 @@ namespace Opm } } + void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + double perm_threshold) + { + const int dim = 3; + const int num_global_cells = grid.cartdims[0]*grid.cartdims[1]*grid.cartdims[2]; + const int nc = grid.number_of_cells; + + assert(num_global_cells > 0); + + permeability_.assign(dim * dim * nc, 0.0); + + std::vector*> tensor; + tensor.reserve(10); + + const std::vector zero(num_global_cells, 0.0); + tensor.push_back(&zero); + + std::array kmap; + PermeabilityKind pkind = fillTensor(newParserDeck, tensor, kmap); + if (pkind == Invalid) { + OPM_THROW(std::runtime_error, "Invalid permeability field."); + } + + // Assign permeability values only if such values are + // given in the input deck represented by 'newParserDeck'. In + // other words: Don't set any (arbitrary) default values. + // It is infinitely better to experience a reproducible + // crash than subtle errors resulting from a (poorly + // chosen) default value... + // + if (tensor.size() > 1) { + const int* gc = grid.global_cell; + int off = 0; + + for (int c = 0; c < nc; ++c, off += dim*dim) { + // SharedPermTensor K(dim, dim, &permeability_[off]); + int kix = 0; + const int glob = (gc == NULL) ? c : gc[c]; + + for (int i = 0; i < dim; ++i) { + for (int j = 0; j < dim; ++j, ++kix) { + // K(i,j) = (*tensor[kmap[kix]])[glob]; + permeability_[off + kix] = (*tensor[kmap[kix]])[glob]; + } + // K(i,i) = std::max(K(i,i), perm_threshold); + permeability_[off + 3*i + i] = std::max(permeability_[off + 3*i + i], perm_threshold); + } + + permfield_valid_[c] = std::vector::value_type(1); + } + } + } namespace { @@ -204,6 +288,72 @@ namespace Opm return retval; } + /// @brief + /// Classify and verify a given permeability specification + /// from a structural point of view. In particular, we + /// verify that there are no off-diagonal permeability + /// components such as @f$k_{xy}@f$ unless the + /// corresponding diagonal components are known as well. + /// + /// @param newParserDeck [in] + /// An Eclipse data parser capable of answering which + /// permeability components are present in a given input + /// deck. + /// + /// @return + /// An enum value with the following possible values: + /// ScalarPerm only one component was given. + /// DiagonalPerm more than one component given. + /// TensorPerm at least one cross-component given. + /// None no components given. + /// Invalid invalid set of components given. + PermeabilityKind classifyPermeability(Opm::DeckConstPtr newParserDeck) + { + const bool xx = newParserDeck->hasKeyword("PERMX" ); + const bool xy = newParserDeck->hasKeyword("PERMXY"); + const bool xz = newParserDeck->hasKeyword("PERMXZ"); + + const bool yx = newParserDeck->hasKeyword("PERMYX"); + const bool yy = newParserDeck->hasKeyword("PERMY" ); + const bool yz = newParserDeck->hasKeyword("PERMYZ"); + + const bool zx = newParserDeck->hasKeyword("PERMZX"); + const bool zy = newParserDeck->hasKeyword("PERMZY"); + const bool zz = newParserDeck->hasKeyword("PERMZ" ); + + int num_cross_comp = xy + xz + yx + yz + zx + zy; + int num_comp = xx + yy + zz + num_cross_comp; + PermeabilityKind retval = None; + if (num_cross_comp > 0) { + retval = TensorPerm; + } else { + if (num_comp == 1) { + retval = ScalarPerm; + } else if (num_comp >= 2) { + retval = DiagonalPerm; + } + } + + bool ok = true; + if (num_comp > 0) { + // At least one tensor component specified on input. + // Verify that any remaining components are OK from a + // structural point of view. In particular, there + // must not be any cross-components (e.g., k_{xy}) + // unless the corresponding diagonal component (e.g., + // k_{xx}) is present as well... + // + ok = xx || !(xy || xz || yx || zx) ; + ok = ok && (yy || !(yx || yz || xy || zy)); + ok = ok && (zz || !(zx || zy || xz || yz)); + } + if (!ok) { + retval = Invalid; + } + + return retval; + } + /// @brief /// Copy isotropic (scalar) permeability to other diagonal @@ -333,6 +483,106 @@ namespace Opm return kind; } + /// @brief + /// Extract pointers to appropriate tensor components from + /// input deck. The permeability tensor is, generally, + /// @code + /// [ kxx kxy kxz ] + /// K = [ kyx kyy kyz ] + /// [ kzx kzy kzz ] + /// @endcode + /// We store these values in a linear array using natural + /// ordering with the column index cycling the most rapidly. + /// In particular we use the representation + /// @code + /// [ 0 1 2 3 4 5 6 7 8 ] + /// K = [ kxx, kxy, kxz, kyx, kyy, kyz, kzx, kzy, kzz ] + /// @endcode + /// Moreover, we explicitly enforce symmetric tensors by + /// assigning + /// @code + /// 3 1 6 2 7 5 + /// kyx = kxy, kzx = kxz, kzy = kyz + /// @endcode + /// However, we make no attempt at enforcing positive + /// definite tensors. + /// + /// @param [in] parser + /// An Eclipse data parser capable of answering which + /// permeability components are present in a given input + /// deck as well as retrieving the numerical value of each + /// permeability component in each grid cell. + /// + /// @param [out] tensor + /// @param [out] kmap + PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck, + std::vector*>& tensor, + std::array& kmap) + { + PermeabilityKind kind = classifyPermeability(newParserDeck); + if (kind == Invalid) { + OPM_THROW(std::runtime_error, "Invalid set of permeability fields given."); + } + assert(tensor.size() == 1); + for (int i = 0; i < 9; ++i) { kmap[i] = 0; } + + enum { xx, xy, xz, // 0, 1, 2 + yx, yy, yz, // 3, 4, 5 + zx, zy, zz }; // 6, 7, 8 + + // ----------------------------------------------------------- + // 1st row: [kxx, kxy, kxz] + if (newParserDeck->hasKeyword("PERMX" )) { + kmap[xx] = tensor.size(); + tensor.push_back(&newParserDeck->getKeyword("PERMX")->getSIDoubleData()); + + setScalarPermIfNeeded(kmap, xx, yy, zz); + } + if (newParserDeck->hasKeyword("PERMXY")) { + kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMXY")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMXZ")) { + kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMXZ")->getSIDoubleData()); + } + + // ----------------------------------------------------------- + // 2nd row: [kyx, kyy, kyz] + if (newParserDeck->hasKeyword("PERMYX")) { + kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMYX")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMY" )) { + kmap[yy] = tensor.size(); + tensor.push_back(&newParserDeck->getKeyword("PERMY")->getSIDoubleData()); + + setScalarPermIfNeeded(kmap, yy, zz, xx); + } + if (newParserDeck->hasKeyword("PERMYZ")) { + kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMYZ")->getSIDoubleData()); + } + + // ----------------------------------------------------------- + // 3rd row: [kzx, kzy, kzz] + if (newParserDeck->hasKeyword("PERMZX")) { + kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMZX")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMZY")) { + kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry. + tensor.push_back(&newParserDeck->getKeyword("PERMZY")->getSIDoubleData()); + } + if (newParserDeck->hasKeyword("PERMZ" )) { + kmap[zz] = tensor.size(); + tensor.push_back(&newParserDeck->getKeyword("PERMZ")->getSIDoubleData()); + + setScalarPermIfNeeded(kmap, zz, xx, yy); + } + return kind; + } + } // anonymous namespace } // namespace Opm diff --git a/opm/core/props/rock/RockFromDeck.hpp b/opm/core/props/rock/RockFromDeck.hpp index e277ee2a..b5e0fb4e 100644 --- a/opm/core/props/rock/RockFromDeck.hpp +++ b/opm/core/props/rock/RockFromDeck.hpp @@ -22,6 +22,9 @@ #include + +#include + #include struct UnstructuredGrid; @@ -43,6 +46,14 @@ namespace Opm void init(const EclipseGridParser& deck, const UnstructuredGrid& grid); + /// Initialize from deck and grid. + /// \param newParserDeck Deck produced by the opm-parser code + /// \param grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + void init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); + /// \return D, the number of spatial dimensions. Always 3 for deck input. int numDimensions() const { @@ -72,9 +83,14 @@ namespace Opm private: void assignPorosity(const EclipseGridParser& parser, const UnstructuredGrid& grid); + void assignPorosity(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); void assignPermeability(const EclipseGridParser& parser, const UnstructuredGrid& grid, const double perm_threshold); + void assignPermeability(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + double perm_threshold); std::vector porosity_; std::vector permeability_; diff --git a/opm/core/props/satfunc/SatFuncBase.hpp b/opm/core/props/satfunc/SatFuncBase.hpp index 47bd31d1..268ecd9d 100644 --- a/opm/core/props/satfunc/SatFuncBase.hpp +++ b/opm/core/props/satfunc/SatFuncBase.hpp @@ -22,6 +22,11 @@ #include #include #include + +#include +#include +#include + #include namespace Opm @@ -78,6 +83,10 @@ namespace Opm const int table_num, const PhaseUsage phase_usg, const int samples); + void init(Opm::DeckConstPtr newParserDeck, + const int table_num, + const PhaseUsage phase_usg, + const int samples); void updateSatHyst(const double* s, const EPSTransforms* epst, const EPSTransforms* epst_hyst, @@ -248,6 +257,138 @@ namespace Opm } } + template + void SatFuncBase::init(Opm::DeckConstPtr newParserDeck, + const int table_num, + const PhaseUsage phase_usg, + const int samples) + { + phase_usage = phase_usg; + double swco = 0.0; + double swmax = 1.0; + if (phase_usage.phase_used[Aqua]) { + Opm::SwofTable swof(newParserDeck->getKeyword("SWOF"), table_num); + const std::vector& sw = swof.getSwColumn(); + const std::vector& krw = swof.getKrwColumn(); + const std::vector& krow = swof.getKrowColumn(); + const std::vector& pcow = swof.getPcowColumn(); + if (krw.front() != 0.0 || krow.back() != 0.0) { + OPM_THROW(std::runtime_error, "Error SWOF data - non-zero krw(swco) and/or krow(1-sor)"); + } + + // Extend the tables with constant values such that the + // derivatives at the endpoints are zero + int n = sw.size(); + std::vector sw_ex(n+2); + std::vector krw_ex(n+2); + std::vector krow_ex(n+2); + std::vector pcow_ex(n+2); + + extendTable(sw,sw_ex,1); + extendTable(krw,krw_ex,0); + extendTable(krow,krow_ex,0); + extendTable(pcow,pcow_ex,0); + + initializeTableType(krw_,sw_ex, krw_ex, samples); + initializeTableType(krow_,sw_ex, krow_ex, samples); + initializeTableType(pcow_,sw_ex, pcow_ex, samples); + + krocw_ = krow[0]; // At connate water -> ecl. SWOF + swco = sw[0]; + smin_[phase_usage.phase_pos[Aqua]] = sw[0]; + swmax = sw.back(); + smax_[phase_usage.phase_pos[Aqua]] = sw.back(); + + krwmax_ = krw.back(); + kromax_ = krow.front(); + swcr_ = swmax; + sowcr_ = 1.0 - swco; + krwr_ = krw.back(); + krorw_ = krow.front(); + for (std::vector::size_type i=1; i 0.0) { + swcr_ = sw[i-1]; + krorw_ = krow[i-1]; + break; + } + } + for (std::vector::size_type i=sw.size()-1; i>=1; --i) { + if (krow[i-1]> 0.0) { + sowcr_ = 1.0 - sw[i]; + krwr_ = krw[i]; + break; + } + } + } + if (phase_usage.phase_used[Vapour]) { + Opm::SgofTable sgof(newParserDeck->getKeyword("SGOF"), table_num); + const std::vector& sg = sgof.getSgColumn(); + const std::vector& krg = sgof.getKrgColumn(); + const std::vector& krog = sgof.getKrogColumn(); + const std::vector& pcog = sgof.getPcogColumn(); + + // Extend the tables with constant values such that the + // derivatives at the endpoints are zero + int n = sg.size(); + std::vector sg_ex(n+2); + std::vector krg_ex(n+2); + std::vector krog_ex(n+2); + std::vector pcog_ex(n+2); + + extendTable(sg,sg_ex,1); + extendTable(krg,krg_ex,0); + extendTable(krog,krog_ex,0); + extendTable(pcog,pcog_ex,0); + + initializeTableType(krg_,sg_ex, krg_ex, samples); + initializeTableType(krog_,sg_ex, krog_ex, samples); + initializeTableType(pcog_,sg_ex, pcog_ex, samples); + + smin_[phase_usage.phase_pos[Vapour]] = sg[0]; + if (std::fabs(sg.back() + swco - 1.0) > 1e-3) { + OPM_THROW(std::runtime_error, "Gas maximum saturation in SGOF table = " << sg.back() << + ", should equal (1.0 - connate water sat) = " << (1.0 - swco)); + } + smax_[phase_usage.phase_pos[Vapour]] = sg.back(); + smin_[phase_usage.phase_pos[Vapour]] = sg.front(); + krgmax_ = krg.back(); + + sgcr_ = sg.front(); + sogcr_ = 1.0 - sg.back(); + krgr_ = krg.back(); + krorg_ = krg.front(); + for (std::vector::size_type i=1; i 0.0) { + sgcr_ = sg[i-1]; + krorg_ = krog[i-1]; + break; + } + } + for (std::vector::size_type i=sg.size()-1; i>=1; --i) { + if (krog[i-1]> 0.0) { + sogcr_ = 1.0 - sg[i]; + krgr_ = krg[i]; + break; + } + } + + } + + if (phase_usage.phase_used[Vapour] && phase_usage.phase_used[Aqua]) { + sowcr_ -= smin_[phase_usage.phase_pos[Vapour]]; + sogcr_ -= smin_[phase_usage.phase_pos[Aqua]]; + smin_[phase_usage.phase_pos[Liquid]] = 0.0; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]] + - smin_[phase_usage.phase_pos[Vapour]]; // First entry in SGOF-table supposed to be zero anyway ... + } else if (phase_usage.phase_used[Aqua]) { + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Aqua]]; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Aqua]]; + } else if (phase_usage.phase_used[Vapour]) { + smin_[phase_usage.phase_pos[Liquid]] = 1.0 - smax_[phase_usage.phase_pos[Vapour]]; + smax_[phase_usage.phase_pos[Liquid]] = 1.0 - smin_[phase_usage.phase_pos[Vapour]]; + } + } + template void SatFuncBase::updateSatHyst(const double* s, const EPSTransforms* epst, diff --git a/opm/core/props/satfunc/SatFuncSimple.cpp b/opm/core/props/satfunc/SatFuncSimple.cpp index 83a377ff..940f12de 100644 --- a/opm/core/props/satfunc/SatFuncSimple.cpp +++ b/opm/core/props/satfunc/SatFuncSimple.cpp @@ -36,7 +36,7 @@ namespace Opm void SatFuncBase >::initializeTableType(NonuniformTableLinear & table, const std::vector& arg, const std::vector& value, - const int /*samples*/) + const int samples) { table = NonuniformTableLinear(arg, value); } @@ -84,7 +84,7 @@ namespace Opm } } - double EPSTransforms::Transform::scaleSatDeriv(double s, double /*s_r*/, double /*s_cr*/, double /*s_max*/) const + double EPSTransforms::Transform::scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const { if (doNotScale) { return 1.0; diff --git a/opm/core/props/satfunc/SatFuncSimple.hpp b/opm/core/props/satfunc/SatFuncSimple.hpp index f2d58b5a..b1781ded 100644 --- a/opm/core/props/satfunc/SatFuncSimple.hpp +++ b/opm/core/props/satfunc/SatFuncSimple.hpp @@ -31,23 +31,23 @@ namespace Opm void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalPc(const double* s, double* pc) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - - void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/) const + + void evalKr(const double* s, double* kr, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} - void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const + void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const; - void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} - void evalPc(const double* /*s*/, double* /*pc*/, const EPSTransforms* /*epst*/) const + void evalPc(const double* s, double* pc, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} - void evalPcDeriv(const double* /*s*/, double* /*pc*/, double* /*dpcds*/, const EPSTransforms* /*epst*/) const + void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");} private: }; - + typedef SatFuncSimple > SatFuncSimpleUniform; typedef SatFuncSimple > SatFuncSimpleNonuniform; @@ -185,7 +185,7 @@ namespace Opm double dkrww = _dsdsw*epst->wat.scaleKrDeriv(s[wpos], this->krw_.derivative(_sw)); double krg = epst->gas.scaleKr(s[gpos], this->krg_(_sg), this->krgr_); double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(s[gpos], this->krg_.derivative(_sg)); - // TODO Check the arguments to the krow- and krog-tables below... + // TODO Check the arguments to the krow- and krog-tables below... double krow = epst->watoil.scaleKr(1.0-s[wpos]-s[gpos], this->krow_(1.0-_sow-this->smin_[gpos]), this->krorw_); // ???? double dkrow = _dsdsow*epst->watoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krow_.derivative(1.0-_sow-this->smin_[gpos])); // ???? //double krog = epst->gasoil.scaleKr(this->krog_(1.0-_sog-this->smin_[wpos]), 1.0-s[wpos]-s[gpos], this->krorg_); // ???? diff --git a/opm/core/props/satfunc/SatFuncStone2.hpp b/opm/core/props/satfunc/SatFuncStone2.hpp index e22f1e9a..4eb5e2fd 100644 --- a/opm/core/props/satfunc/SatFuncStone2.hpp +++ b/opm/core/props/satfunc/SatFuncStone2.hpp @@ -31,25 +31,25 @@ namespace Opm void evalKrDeriv(const double* s, double* kr, double* dkrds) const; void evalPc(const double* s, double* pc) const; void evalPcDeriv(const double* s, double* pc, double* dpcds) const; - - void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/) const + + void evalKr(const double* s, double* kr, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} - void evalKr(const double* /*s*/, double* /*kr*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const + void evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} - void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/) const + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} - void evalKrDeriv(const double* /*s*/, double* /*kr*/, double* /*dkrds*/, const EPSTransforms* /*epst*/, const EPSTransforms* /*epst_hyst*/, const SatHyst* /*sat_hyst*/) const + void evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} - void evalPc(const double* /*s*/, double* /*pc*/, const EPSTransforms* /*epst*/) const + void evalPc(const double* s, double* pc, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} - void evalPcDeriv(const double* /*s*/, double* /*pc*/, double* /*dpcds*/, const EPSTransforms* /*epst*/) const + void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const {OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");} private: }; - + typedef SatFuncStone2 > SatFuncStone2Uniform; typedef SatFuncStone2 > SatFuncStone2Nonuniform; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp index c321a42a..2c570ffd 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck.hpp @@ -27,6 +27,9 @@ #include #include #include + +#include + #include struct UnstructuredGrid; @@ -60,6 +63,17 @@ namespace Opm const UnstructuredGrid& grid, const int samples); + /// Initialize from deck and grid. + /// \param[in] deck Deck input parser + /// \param[in] grid Grid to which property object applies, needed for the + /// mapping from cell indices (typically from a processed grid) + /// to logical cartesian indices consistent with the deck. + /// \param[in] samples Number of uniform sample points for saturation tables. + /// NOTE: samples will only be used with the SatFuncSetUniform template argument. + void init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const int samples); + /// \return P, the number of phases. int numPhases() const; @@ -132,6 +146,14 @@ namespace Opm const UnstructuredGrid& grid, const std::string& keyword, std::vector& scaleparam); + void initEPS(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); + void initEPSHyst(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid); + void initEPSKey(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam); void initEPSParam(const int cell, EPSTransforms::Transform& data, const bool oil, @@ -149,6 +171,11 @@ namespace Opm const std::vector& s0, const std::vector& krsr, const std::vector& krmax); + + bool columnIsMasked_(Opm::DeckConstPtr newParserDeck, + const std::string &keywordName, + int columnIdx) + { return newParserDeck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; } }; diff --git a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp index c3006993..d9a3f32e 100644 --- a/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp +++ b/opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp @@ -26,6 +26,11 @@ #include #include +#include +#include +#include +#include + #include namespace Opm @@ -140,6 +145,123 @@ namespace Opm } + /// Initialize from deck. + template + void SaturationPropsFromDeck::init(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const int samples) + { + phase_usage_ = phaseUsageFromDeck(newParserDeck); + + // Extract input data. + // Oil phase should be active. + if (!phase_usage_.phase_used[Liquid]) { + OPM_THROW(std::runtime_error, "SaturationPropsFromDeck::init() -- oil phase must be active."); + } + + // Obtain SATNUM, if it exists, and create cell_to_func_. + // Otherwise, let the cell_to_func_ mapping be just empty. + int satfuncs_expected = 1; + if (newParserDeck->hasKeyword("SATNUM")) { + const std::vector& satnum = newParserDeck->getKeyword("SATNUM")->getIntData(); + satfuncs_expected = *std::max_element(satnum.begin(), satnum.end()); + const int num_cells = grid.number_of_cells; + cell_to_func_.resize(num_cells); + const int* gc = grid.global_cell; + for (int cell = 0; cell < num_cells; ++cell) { + const int newParserDeck_pos = (gc == NULL) ? cell : gc[cell]; + cell_to_func_[cell] = satnum[newParserDeck_pos] - 1; + } + } + + // Find number of tables, check for consistency. + enum { Uninitialized = -1 }; + int num_tables = Uninitialized; + if (phase_usage_.phase_used[Aqua]) { + num_tables = newParserDeck->getKeyword("SWOF")->size(); + if (num_tables < satfuncs_expected) { + OPM_THROW(std::runtime_error, "Found " << num_tables << " SWOF tables, SATNUM specifies at least " << satfuncs_expected); + } + } + if (phase_usage_.phase_used[Vapour]) { + int num_sgof_tables = newParserDeck->getKeyword("SGOF")->size(); + if (num_sgof_tables < satfuncs_expected) { + OPM_THROW(std::runtime_error, "Found " << num_tables << " SGOF tables, SATNUM specifies at least " << satfuncs_expected); + } + if (num_tables == Uninitialized) { + num_tables = num_sgof_tables; + } else if (num_tables != num_sgof_tables) { + OPM_THROW(std::runtime_error, "Inconsistent number of tables in SWOF and SGOF."); + } + } + + // Initialize tables. + satfuncset_.resize(num_tables); + for (int table = 0; table < num_tables; ++table) { + satfuncset_[table].init(newParserDeck, table, phase_usage_, samples); + } + + // Saturation table scaling + do_hyst_ = false; + do_eps_ = false; + do_3pt_ = false; + if (newParserDeck->hasKeyword("ENDSCALE")) { + Opm::EndscaleWrapper endscale(newParserDeck->getKeyword("ENDSCALE")); + if (endscale.directionSwitch() != std::string("NODIR")) { + OPM_THROW(std::runtime_error, + "SaturationPropsFromDeck::init() -- ENDSCALE: " + "Currently only 'NODIR' accepted."); + } + if (endscale.isReversible()) { + OPM_THROW(std::runtime_error, + "SaturationPropsFromDeck::init() -- ENDSCALE: " + "Currently only 'REVERS' accepted."); + } + if (newParserDeck->hasKeyword("SCALECRS")) { + Opm::ScalecrsWrapper scalecrs(newParserDeck->getKeyword("SCALECRS")); + if (scalecrs.isEnabled()) { + do_3pt_ = true; + } + } + do_eps_ = true; + + initEPS(newParserDeck, grid); + + // For now, a primitive detection of hysteresis. TODO: SATOPTS HYSTER/ and EHYSTR + do_hyst_ = + newParserDeck->hasKeyword("ISWL") + || newParserDeck->hasKeyword("ISWU") + || newParserDeck->hasKeyword("ISWCR") + || newParserDeck->hasKeyword("ISGL") + || newParserDeck->hasKeyword("ISGU") + || newParserDeck->hasKeyword("ISGCR") + || newParserDeck->hasKeyword("ISOWCR") + || newParserDeck->hasKeyword("ISOGCR"); + if (do_hyst_) { + if (newParserDeck->hasKeyword("KRW") + || newParserDeck->hasKeyword("KRG") + || newParserDeck->hasKeyword("KRO") + || newParserDeck->hasKeyword("KRWR") + || newParserDeck->hasKeyword("KRGR") + || newParserDeck->hasKeyword("KRORW") + || newParserDeck->hasKeyword("KRORG") + || newParserDeck->hasKeyword("IKRW") + || newParserDeck->hasKeyword("IKRG") + || newParserDeck->hasKeyword("IKRO") + || newParserDeck->hasKeyword("IKRWR") + || newParserDeck->hasKeyword("IKRGR") + || newParserDeck->hasKeyword("IKRORW") + || newParserDeck->hasKeyword("IKRORG") ) { + OPM_THROW(std::runtime_error, + "SaturationPropsFromDeck::init() -- ENDSCALE: " + "Currently hysteresis and relperm value scaling " + "cannot be combined."); + } + initEPSHyst(newParserDeck, grid); + } + } + } + /// \return P, the number of phases. @@ -383,6 +505,69 @@ namespace Opm } } + // Initialize saturation scaling parameters + template + void SaturationPropsFromDeck::initEPS(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + std::vector swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr; + std::vector krw, krg, kro, krwr, krgr, krorw, krorg; + // Initialize saturation scaling parameter + initEPSKey(newParserDeck, grid, std::string("SWL"), swl); + initEPSKey(newParserDeck, grid, std::string("SWU"), swu); + initEPSKey(newParserDeck, grid, std::string("SWCR"), swcr); + initEPSKey(newParserDeck, grid, std::string("SGL"), sgl); + initEPSKey(newParserDeck, grid, std::string("SGU"), sgu); + initEPSKey(newParserDeck, grid, std::string("SGCR"), sgcr); + initEPSKey(newParserDeck, grid, std::string("SOWCR"), sowcr); + initEPSKey(newParserDeck, grid, std::string("SOGCR"), sogcr); + initEPSKey(newParserDeck, grid, std::string("KRW"), krw); + initEPSKey(newParserDeck, grid, std::string("KRG"), krg); + initEPSKey(newParserDeck, grid, std::string("KRO"), kro); + initEPSKey(newParserDeck, grid, std::string("KRWR"), krwr); + initEPSKey(newParserDeck, grid, std::string("KRGR"), krgr); + initEPSKey(newParserDeck, grid, std::string("KRORW"), krorw); + initEPSKey(newParserDeck, grid, std::string("KRORG"), krorg); + + eps_transf_.resize(grid.number_of_cells); + + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; + const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + if (oilWater) { + // ### krw + initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw); + // ### krow + initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro); + } else if (oilGas) { + // ### krg + initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg); + // ### krog + initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro); + } else if (threephase) { + // ### krw + initEPSParam(cell, eps_transf_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, swl, swcr, swu, sowcr, sgl, krwr, krw); + // ### krow + initEPSParam(cell, eps_transf_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, swl, sowcr, swl, swcr, sgl, krorw, kro); + // ### krg + initEPSParam(cell, eps_transf_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, sgl, sgcr, sgu, sogcr, swl, krgr, krg); + // ### krog + initEPSParam(cell, eps_transf_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, sgl, sogcr, sgl, sgcr, swl, krorg, kro); + } + } + } // Initialize hysteresis saturation scaling parameters template @@ -450,6 +635,71 @@ namespace Opm } } + // Initialize hysteresis saturation scaling parameters + template + void SaturationPropsFromDeck::initEPSHyst(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid) + { + std::vector iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr; + std::vector ikrw, ikrg, ikro, ikrwr, ikrgr, ikrorw, ikrorg; + // Initialize hysteresis saturation scaling parameters + initEPSKey(newParserDeck, grid, std::string("ISWL"), iswl); + initEPSKey(newParserDeck, grid, std::string("ISWU"), iswu); + initEPSKey(newParserDeck, grid, std::string("ISWCR"), iswcr); + initEPSKey(newParserDeck, grid, std::string("ISGL"), isgl); + initEPSKey(newParserDeck, grid, std::string("ISGU"), isgu); + initEPSKey(newParserDeck, grid, std::string("ISGCR"), isgcr); + initEPSKey(newParserDeck, grid, std::string("ISOWCR"), isowcr); + initEPSKey(newParserDeck, grid, std::string("ISOGCR"), isogcr); + initEPSKey(newParserDeck, grid, std::string("IKRW"), ikrw); + initEPSKey(newParserDeck, grid, std::string("IKRG"), ikrg); + initEPSKey(newParserDeck, grid, std::string("IKRO"), ikro); + initEPSKey(newParserDeck, grid, std::string("IKRWR"), ikrwr); + initEPSKey(newParserDeck, grid, std::string("IKRGR"), ikrgr); + initEPSKey(newParserDeck, grid, std::string("IKRORW"), ikrorw); + initEPSKey(newParserDeck, grid, std::string("IKRORG"), ikrorg); + + eps_transf_hyst_.resize(grid.number_of_cells); + sat_hyst_.resize(grid.number_of_cells); + + const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + const int gpos = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + const bool oilWater = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && !phase_usage_.phase_used[Vapour]; + const bool oilGas = !phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + const bool threephase = phase_usage_.phase_used[Aqua] && phase_usage_.phase_used[Liquid] && phase_usage_.phase_used[Vapour]; + + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + if (oilWater) { + // ### krw + initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], + funcForCell(cell).sowcr_, -1.0, funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw); + // ### krow + initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], + funcForCell(cell).swcr_, -1.0, funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro); + } else if (oilGas) { + // ### krg + initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], + funcForCell(cell).sogcr_, -1.0, funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg); + // ### krog + initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], + funcForCell(cell).sgcr_, -1.0, funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro); + } else if (threephase) { + // ### krw + initEPSParam(cell, eps_transf_hyst_[cell].wat, false, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, funcForCell(cell).smax_[wpos], funcForCell(cell).sowcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krwr_, funcForCell(cell).krwmax_, iswl, iswcr, iswu, isowcr, isgl, ikrwr, ikrw); + // ### krow + initEPSParam(cell, eps_transf_hyst_[cell].watoil, true, 0.0, funcForCell(cell).sowcr_, funcForCell(cell).smin_[wpos], funcForCell(cell).swcr_, + funcForCell(cell).smin_[gpos], funcForCell(cell).krorw_, funcForCell(cell).kromax_, iswl, isowcr, iswl, iswcr, isgl, ikrorw, ikro); + // ### krg + initEPSParam(cell, eps_transf_hyst_[cell].gas, false, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, funcForCell(cell).smax_[gpos], funcForCell(cell).sogcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krgr_, funcForCell(cell).krgmax_, isgl, isgcr, isgu, isogcr, iswl, ikrgr, ikrg); + // ### krog + initEPSParam(cell, eps_transf_hyst_[cell].gasoil, true, 0.0, funcForCell(cell).sogcr_, funcForCell(cell).smin_[gpos], funcForCell(cell).sgcr_, + funcForCell(cell).smin_[wpos], funcForCell(cell).krorg_, funcForCell(cell).kromax_, isgl, isogcr, isgl, isgcr, iswl, ikrorg, ikro); + } + } + } + // Initialize saturation scaling parameter template void SaturationPropsFromDeck::initEPSKey(const EclipseGridParser& deck, @@ -624,6 +874,213 @@ namespace Opm } } + // Initialize saturation scaling parameter + template + void SaturationPropsFromDeck::initEPSKey(Opm::DeckConstPtr newParserDeck, + const UnstructuredGrid& grid, + const std::string& keyword, + std::vector& scaleparam) + { + const bool useAqua = phase_usage_.phase_used[Aqua]; + const bool useLiquid = phase_usage_.phase_used[Liquid]; + const bool useVapour = phase_usage_.phase_used[Vapour]; + bool useKeyword = newParserDeck->hasKeyword(keyword); + bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD"); + bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD"); + int itab = 0; + std::vector > > table_dummy; + std::vector > >& table = table_dummy; + + // Active keyword assigned default values for each cell (in case of possible box-wise assignment) + int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua]; + int phase_pos_vapour = phase_usage_.phase_pos[BlackoilPhases::Vapour]; + if ((keyword[0] == 'S' && (useKeyword || hasENPTVD)) || (keyword[1] == 'S' && useKeyword) ) { + if (keyword == std::string("SWL") || keyword == std::string("ISWL") ) { + if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 0))) { + itab = 1; + scaleparam.resize(grid.number_of_cells); + for (int i=0; i 0) { + Opm::EnptvdTable enptvd(newParserDeck->getKeyword("ENPTVD")); + table.resize(1); // only one region supported so far + for (unsigned i = 0; i < table.size(); ++i) { + table[i].resize(9); + + for (int k = 0; k < enptvd.numRows(); ++k) { + table[i][0] = enptvd.getDepthColumn(); + table[i][1] = enptvd.getSwcoColumn(); + table[i][2] = enptvd.getSwcritColumn(); + table[i][3] = enptvd.getSwmaxColumn(); + table[i][4] = enptvd.getSgcoColumn(); + table[i][5] = enptvd.getSgcritColumn(); + table[i][6] = enptvd.getSgmaxColumn(); + table[i][7] = enptvd.getSowcritColumn(); + table[i][8] = enptvd.getSogcritColumn(); + } + } + } + } else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) { + if (keyword == std::string("KRW") || keyword == std::string("IKRW") ) { + if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 0))) { + itab = 1; + scaleparam.resize(grid.number_of_cells); + for (int i=0; i 0) { + Opm::EnkrvdTable enkrvd(newParserDeck->getKeyword("ENKRVD")); + table.resize(1); // only one region supported so far + for (unsigned i = 0; i < table.size(); ++i) { + table[i].resize(8); + + for (int k = 0; k < enkrvd.numRows(); ++k) { + table[i][0] = enkrvd.getDepthColumn(); + table[i][1] = enkrvd.getKrwmaxColumn(); + table[i][2] = enkrvd.getKrgmaxColumn(); + table[i][3] = enkrvd.getKromaxColumn(); + table[i][4] = enkrvd.getKrwcritColumn(); + table[i][5] = enkrvd.getKrgcritColumn(); + table[i][6] = enkrvd.getKrocritgColumn(); + table[i][7] = enkrvd.getKrocritwColumn(); + } + } + } + } + + if (scaleparam.empty()) { + return; + } else if (useKeyword) { + // Keyword values from deck + std::cout << "--- Scaling parameter '" << keyword << "' assigned." << std::endl; + const int* gc = grid.global_cell; + const std::vector& val = newParserDeck->getKeyword(keyword)->getSIDoubleData(); + for (int c = 0; c < int(scaleparam.size()); ++c) { + const int deck_pos = (gc == NULL) ? c : gc[c]; + scaleparam[c] = val[deck_pos]; + } + } else { + // TODO for new parser + /* + std::cout << "--- Scaling parameter '" << keyword << "' assigned via "; + if (keyword[0] == 'S') + newParserDeck.getENPTVD().write(std::cout); + else + newParserDeck.getENKRVD().write(std::cout); + */ + const double* cc = grid.cell_centroids; + const int dim = grid.dimensions; + for (int cell = 0; cell < grid.number_of_cells; ++cell) { + int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell]; + if (table[itab][jtab][0] != -1.0) { + std::vector& depth = table[0][jtab]; + std::vector& val = table[itab][jtab]; + double zc = cc[dim*cell+dim-1]; + if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval + scaleparam[cell] = linearInterpolation(depth, val, zc); + } + } + } + } + } + // Saturation scaling template void SaturationPropsFromDeck::initEPSParam(const int cell, diff --git a/opm/core/simulator/SimulatorTimer.cpp b/opm/core/simulator/SimulatorTimer.cpp index f3bef455..fddc3aa8 100644 --- a/opm/core/simulator/SimulatorTimer.cpp +++ b/opm/core/simulator/SimulatorTimer.cpp @@ -57,10 +57,21 @@ namespace Opm start_date_ = deck.getStartDate(); } + /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap + void SimulatorTimer::init(Opm::TimeMapConstPtr timeMap, + int timeStepIdx) + { + timeMap_ = timeMap; + current_step_ = timeStepIdx; + } + /// Total number of steps. int SimulatorTimer::numSteps() const { - return timesteps_.size(); + if (timeMap_) + return timeMap_->numTimesteps(); + else + return timesteps_.size(); } /// Current step number. @@ -72,13 +83,14 @@ namespace Opm /// Set current step number. void SimulatorTimer::setCurrentStepNum(int step) { - if (current_step_ < 0 || current_step_ > int(timesteps_.size())) { + if (current_step_ < 0 || current_step_ > int(numSteps())) { // Note that we do allow current_step_ == timesteps_.size(), // that is the done() state. OPM_THROW(std::runtime_error, "Trying to set invalid step number: " << step); } current_step_ = step; - current_time_ = std::accumulate(timesteps_.begin(), timesteps_.begin() + step, 0.0); + if (timeMap_) + current_time_ = std::accumulate(timesteps_.begin(), timesteps_.begin() + step, 0.0); } @@ -86,25 +98,37 @@ namespace Opm double SimulatorTimer::currentStepLength() const { assert(!done()); - return timesteps_[current_step_]; + if (timeMap_) + return timeMap_->getTimeStepLength(current_step_); + else + return timesteps_[current_step_]; } double SimulatorTimer::stepLengthTaken() const { assert(current_step_ > 0); - return timesteps_[current_step_ - 1]; + if (timeMap_) + return timeMap_->getTimeStepLength(current_step_ - 1); + else + return timesteps_[current_step_ - 1]; } /// Current time. double SimulatorTimer::currentTime() const { - return current_time_; + if (timeMap_) + return timeMap_->getTimePassedUntil(current_step_); + else + return current_time_; } boost::posix_time::ptime SimulatorTimer::currentDateTime() const { - return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ ); + if (timeMap_) + return timeMap_->getStartTime(current_step_); + else + return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ ); } @@ -112,7 +136,10 @@ namespace Opm /// Total time. double SimulatorTimer::totalTime() const { - return total_time_; + if (timeMap_) + return timeMap_->getTotalTime(); + else + return total_time_; } /// Set total time. @@ -121,7 +148,13 @@ namespace Opm /// access to later timesteps. void SimulatorTimer::setTotalTime(double time) { - total_time_ = time; + if (timeMap_) { + // well, what can we do if we use opm-parser's TimeMap? + OPM_THROW(std::logic_error, + "Not implemented: SimulatorTimer::setTotalTime() if using a TimeMap."); + } + else + total_time_ = time; } /// Print a report with current and total time etc. @@ -138,7 +171,8 @@ namespace Opm SimulatorTimer& SimulatorTimer::operator++() { assert(!done()); - current_time_ += timesteps_[current_step_]; + if (!timeMap_) + current_time_ += timesteps_[current_step_]; ++current_step_; return *this; } @@ -146,7 +180,10 @@ namespace Opm /// Return true if op++() has been called numSteps() times. bool SimulatorTimer::done() const { - return int(timesteps_.size()) == current_step_; + if (timeMap_) + return int(timeMap_->numTimesteps()) <= current_step_; + else + return int(timesteps_.size()) == current_step_; } diff --git a/opm/core/simulator/SimulatorTimer.hpp b/opm/core/simulator/SimulatorTimer.hpp index a1d8a270..949859e5 100644 --- a/opm/core/simulator/SimulatorTimer.hpp +++ b/opm/core/simulator/SimulatorTimer.hpp @@ -20,6 +20,8 @@ #ifndef OPM_SIMULATORTIMER_HEADER_INCLUDED #define OPM_SIMULATORTIMER_HEADER_INCLUDED +#include + #include #include #include @@ -47,6 +49,10 @@ namespace Opm /// Note that DATES are folded into TSTEP by the parser. void init(const EclipseGridParser& deck); + /// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap + void init(TimeMapConstPtr timeMap, + int timeStepIdx = 0); + /// Total number of steps. int numSteps() const; @@ -99,6 +105,7 @@ namespace Opm bool done() const; private: + Opm::TimeMapConstPtr timeMap_; std::vector timesteps_; int current_step_; double current_time_; diff --git a/opm/core/simulator/initState_impl.hpp b/opm/core/simulator/initState_impl.hpp index 12dbe7b5..18372573 100644 --- a/opm/core/simulator/initState_impl.hpp +++ b/opm/core/simulator/initState_impl.hpp @@ -32,6 +32,8 @@ #include #include +#include + #include #include @@ -471,6 +473,103 @@ namespace Opm } + /// Initialize a state from input deck. + template + void initStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + Opm::DeckConstPtr newParserDeck, + const double gravity, + State& state) + { + const int num_phases = props.numPhases(); + const PhaseUsage pu = phaseUsageFromDeck(newParserDeck); + if (num_phases != pu.num_phases) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, " + "found " << pu.num_phases << " phases in deck."); + } + state.init(grid, num_phases); + if (newParserDeck->hasKeyword("EQUIL")) { + if (num_phases != 2) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); + } + if (pu.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas)."); + } + // Set saturations depending on oil-water contact. + EquilWrapper equil(newParserDeck->getKeyword("EQUIL")); + if (equil.numRegions() != 1) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet."); + } + const double woc = equil.waterOilContactDepth(0); + initWaterOilContact(grid, props, woc, WaterBelow, state); + // Set pressure depending on densities and depths. + const double datum_z = equil.datumDepth(0); + const double datum_p = equil.datumDepthPressure(0); + initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); + } else if (newParserDeck->hasKeyword("PRESSURE")) { + // Set saturations from SWAT/SGAS, pressure from PRESSURE. + std::vector& s = state.saturation(); + std::vector& p = state.pressure(); + const std::vector& p_deck = newParserDeck->getKeyword("PRESSURE")->getSIDoubleData(); + const int num_cells = grid.number_of_cells; + if (num_phases == 2) { + if (!pu.phase_used[BlackoilPhases::Aqua]) { + // oil-gas: we require SGAS + if (!newParserDeck->hasKeyword("SGAS")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init"); + } + const std::vector& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData(); + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + gpos] = sg_deck[c_deck]; + s[2*c + opos] = 1.0 - sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + // water-oil or water-gas: we require SWAT + if (!newParserDeck->hasKeyword("SWAT")) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init"); + } + const std::vector& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData(); + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int nwpos = (wpos + 1) % 2; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + wpos] = sw_deck[c_deck]; + s[2*c + nwpos] = 1.0 - sw_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } + } else if (num_phases == 3) { + const bool has_swat_sgas = newParserDeck->hasKeyword("SWAT") && newParserDeck->hasKeyword("SGAS"); + if (!has_swat_sgas) { + OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init."); + } + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData(); + const std::vector& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData(); + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[3*c + wpos] = sw_deck[c_deck]; + s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); + s[3*c + gpos] = sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); + } + } else { + OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS."); + } + + // Finally, init face pressures. + initFacePressure(grid, state); + } + /// Initialize a state from input deck. template void initStateFromDeck(const UnstructuredGrid& grid, @@ -568,10 +667,6 @@ namespace Opm initFacePressure(grid, state); } - - - - /// Initialize surface volume from pressure and saturation by z = As. /// Here saturation is used as an intial guess for z in the /// computation of A. @@ -770,6 +865,39 @@ namespace Opm } } + /// Initialize a blackoil state from input deck. + template + void initBlackoilStateFromDeck(const UnstructuredGrid& grid, + const Props& props, + Opm::DeckConstPtr newParserDeck, + const double gravity, + State& state) + { + initStateFromDeck(grid, props, newParserDeck, gravity, state); + if (newParserDeck->hasKeyword("RS")) { + const std::vector& rs_deck = newParserDeck->getKeyword("RS")->getSIDoubleData(); + const int num_cells = grid.number_of_cells; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + state.gasoilratio()[c] = rs_deck[c_deck]; + } + initBlackoilSurfvolUsingRSorRV(grid, props, state); + computeSaturation(props,state); + } else if (newParserDeck->hasKeyword("RV")){ + const std::vector& rv_deck = newParserDeck->getKeyword("RV")->getSIDoubleData(); + const int num_cells = grid.number_of_cells; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + state.rv()[c] = rv_deck[c_deck]; + } + initBlackoilSurfvolUsingRSorRV(grid, props, state); + computeSaturation(props,state); + } + else { + OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field."); + } + } + } // namespace Opm diff --git a/opm/core/well_controls.h b/opm/core/well_controls.h index 70ea1b49..435c7a79 100644 --- a/opm/core/well_controls.h +++ b/opm/core/well_controls.h @@ -34,10 +34,10 @@ enum WellControlType { SURFACE_RATE /**< Well constrained by surface volume flow rate */ }; -struct WellControls; +struct WellControls; bool -well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2); +well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2 , bool verbose); struct WellControls * well_controls_create(void); diff --git a/opm/core/wells.h b/opm/core/wells.h index dfc9be17..05b6ce3e 100644 --- a/opm/core/wells.h +++ b/opm/core/wells.h @@ -263,7 +263,7 @@ struct Wells * clone_wells(const struct Wells *W); bool -wells_equal(const struct Wells *W1, const struct Wells *W2); +wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose); diff --git a/opm/core/wells/WellCollection.cpp b/opm/core/wells/WellCollection.cpp index 350d2ffb..9218db06 100644 --- a/opm/core/wells/WellCollection.cpp +++ b/opm/core/wells/WellCollection.cpp @@ -19,10 +19,66 @@ #include "config.h" #include + +#include +#include + +#include + #include namespace Opm { + void WellCollection::addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage) { + WellsGroupInterface* fieldNode = findNode(fieldGroup->name()); + if (fieldNode) { + OPM_THROW(std::runtime_error, "Trying to add FIELD node, but this already exists. Can only have one FIELD node."); + } + + roots_.push_back(createGroupWellsGroup(fieldGroup, timeStep, phaseUsage)); + } + + void WellCollection::addGroup(GroupConstPtr groupChild, std::string parent_name, + size_t timeStep, const PhaseUsage& phaseUsage) { + WellsGroupInterface* parent = findNode(parent_name); + if (!parent) { + OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent_name << ", but this does not exist in the WellCollection."); + } + + if (findNode(groupChild->name())) { + OPM_THROW(std::runtime_error, "Trying to add child group named " << groupChild->name() << ", but this group is already in the WellCollection."); + + } + + std::shared_ptr child = createGroupWellsGroup(groupChild, timeStep, phaseUsage); + + WellsGroup* parent_as_group = static_cast (parent); + if (!parent_as_group) { + OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent->name() << ", but it's not a group."); + } + parent_as_group->addChild(child); + child->setParent(parent); + } + + void WellCollection::addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage) { + WellsGroupInterface* parent = findNode(wellChild->getGroupName(timeStep)); + if (!parent) { + OPM_THROW(std::runtime_error, "Trying to add well " << wellChild->name() << " Step: " << boost::lexical_cast(timeStep) << " to group named " << wellChild->getGroupName(timeStep) << ", but this group does not exist in the WellCollection."); + } + + std::shared_ptr child = createWellWellsGroup(wellChild, timeStep, phaseUsage); + + WellsGroup* parent_as_group = static_cast (parent); + if (!parent_as_group) { + OPM_THROW(std::runtime_error, "Trying to add well to group named " << wellChild->getGroupName(timeStep) << ", but it's not a group."); + } + parent_as_group->addChild(child); + + leaf_nodes_.push_back(static_cast(child.get())); + + child->setParent(parent); + } + void WellCollection::addChild(const std::string& child_name, const std::string& parent_name, @@ -33,6 +89,7 @@ namespace Opm roots_.push_back(createWellsGroup(parent_name, deck)); parent = roots_[roots_.size() - 1].get(); } + std::shared_ptr child; for (size_t i = 0; i < roots_.size(); ++i) { @@ -47,6 +104,7 @@ namespace Opm break; } } + if (!child.get()) { child = createWellsGroup(child_name, deck); } @@ -65,7 +123,6 @@ namespace Opm } - const std::vector& WellCollection::getLeafNodes() const { return leaf_nodes_; } diff --git a/opm/core/wells/WellCollection.hpp b/opm/core/wells/WellCollection.hpp index 9cdb9e7f..f084a04d 100644 --- a/opm/core/wells/WellCollection.hpp +++ b/opm/core/wells/WellCollection.hpp @@ -23,10 +23,15 @@ #define OPM_WELLCOLLECTION_HPP #include +#include + #include #include #include -#include +#include + +#include +#include namespace Opm { @@ -34,6 +39,14 @@ namespace Opm class WellCollection { public: + + void addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage); + + void addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage); + + void addGroup(GroupConstPtr groupChild, std::string parent_name, + size_t timeStep, const PhaseUsage& phaseUsage); + /// Adds and creates if necessary the child to the collection /// and appends it to parent's children. Also adds and creates the parent /// if necessary. diff --git a/opm/core/wells/WellsGroup.cpp b/opm/core/wells/WellsGroup.cpp index 6c1059a2..3316ccdb 100644 --- a/opm/core/wells/WellsGroup.cpp +++ b/opm/core/wells/WellsGroup.cpp @@ -74,7 +74,7 @@ namespace Opm { return parent_; } - const std::string& WellsGroupInterface::name() + const std::string& WellsGroupInterface::name() const { return name_; } @@ -1141,4 +1141,54 @@ namespace Opm return return_value; } + + std::shared_ptr createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage ) + { + InjectionSpecification injection_specification; + ProductionSpecification production_specification; + if (group->isInjectionGroup(timeStep)) { + injection_specification.injector_type_ = toInjectorType(Phase::PhaseEnum2String(group->getInjectionPhase(timeStep))); + injection_specification.control_mode_ = toInjectionControlMode(GroupInjection::ControlEnum2String(group->getInjectionControlMode(timeStep))); + injection_specification.surface_flow_max_rate_ = group->getSurfaceMaxRate(timeStep); + injection_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep); + injection_specification.reinjection_fraction_target_ = group->getTargetReinjectFraction(timeStep); + injection_specification.voidage_replacment_fraction_ = group->getTargetVoidReplacementFraction(timeStep); + } + else if (group->isProductionGroup(timeStep)) { + production_specification.oil_max_rate_ = group->getOilTargetRate(timeStep); + production_specification.control_mode_ = toProductionControlMode(GroupProduction::ControlEnum2String(group->getProductionControlMode(timeStep))); + production_specification.water_max_rate_ = group->getWaterTargetRate(timeStep); + production_specification.gas_max_rate_ = group->getGasTargetRate(timeStep); + production_specification.liquid_max_rate_ = group->getLiquidTargetRate(timeStep); + production_specification.procedure_ = toProductionProcedure(GroupProductionExceedLimit::ActionEnum2String(group->getProductionExceedLimitAction(timeStep))); + production_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep); + } + + std::shared_ptr wells_group(new WellsGroup(group->name(), production_specification, injection_specification, phase_usage)); + return wells_group; + } + + std::shared_ptr createWellWellsGroup(WellConstPtr well, size_t timeStep, const PhaseUsage& phase_usage ) + { + InjectionSpecification injection_specification; + ProductionSpecification production_specification; + if (well->isInjector(timeStep)) { + injection_specification.BHP_limit_ = well->getBHPLimit(timeStep); + injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(well->getInjectorType(timeStep))); + injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(well->getInjectorControlMode(timeStep))); + injection_specification.surface_flow_max_rate_ = well->getSurfaceInjectionRate(timeStep); + injection_specification.reservoir_flow_max_rate_ = well->getReservoirInjectionRate(timeStep); + production_specification.guide_rate_ = 0.0; // We know we're not a producer + } + else if (well->isProducer(timeStep)) { + production_specification.BHP_limit_ = well->getBHPLimit(timeStep); + production_specification.reservoir_flow_max_rate_ = well->getResVRate(timeStep); + production_specification.oil_max_rate_ = well->getOilRate(timeStep); + production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(well->getProducerControlMode(timeStep))); + production_specification.water_max_rate_ = well->getWaterRate(timeStep); + injection_specification.guide_rate_ = 0.0; // we know we're not an injector + } + std::shared_ptr wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage)); + return wells_group; + } } diff --git a/opm/core/wells/WellsGroup.hpp b/opm/core/wells/WellsGroup.hpp index 143da935..03a61d5e 100644 --- a/opm/core/wells/WellsGroup.hpp +++ b/opm/core/wells/WellsGroup.hpp @@ -25,6 +25,10 @@ #include #include #include + +#include +#include + #include #include @@ -58,7 +62,7 @@ namespace Opm virtual ~WellsGroupInterface(); /// The unique identifier for the well or well group. - const std::string& name(); + const std::string& name() const; /// Production specifications for the well or well group. const ProductionSpecification& prodSpec() const; @@ -403,9 +407,21 @@ namespace Opm /// \param[in] name the name of the wells group. /// \param[in] deck the deck from which to fetch information. std::shared_ptr createWellsGroup(const std::string& name, - const EclipseGridParser& deck); + const EclipseGridParser& deck); + /// Creates the WellsGroupInterface for the given well + /// \param[in] well the Well to construct object for + /// \param[in] timeStep the time step in question + /// \param[in] the phase usage + std::shared_ptr createWellWellsGroup(WellConstPtr well, size_t timeStep, + const PhaseUsage& phase_usage ); + /// Creates the WellsGroupInterface for the given Group + /// \param[in] group the Group to construct object for + /// \param[in] timeStep the time step in question + /// \param[in] the phase usage + std::shared_ptr createGroupWellsGroup(GroupConstPtr group, size_t timeStep, + const PhaseUsage& phase_usage ); } #endif /* OPM_WELLSGROUP_HPP */ diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 007c8b52..64091120 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -18,6 +18,8 @@ */ #include "config.h" + + #include #include #include @@ -28,6 +30,8 @@ #include #include +#include + #include #include #include @@ -43,23 +47,8 @@ namespace { - struct WellData - { - WellType type; - // WellControlType control; - // double target; - double reference_bhp_depth; - // Opm::InjectionSpecification::InjectorType injected_phase; - int welspecsline; - }; - struct PerfData - { - int cell; - double well_index; - }; - namespace ProductionControl { enum Mode { ORAT, WRAT, GRAT, @@ -101,6 +90,33 @@ namespace << control << " in input file"); } } + + + Mode mode(Opm::WellProducer::ControlModeEnum controlMode) + { + switch( controlMode ) { + case Opm::WellProducer::ORAT: + return ORAT; + case Opm::WellProducer::WRAT: + return WRAT; + case Opm::WellProducer::GRAT: + return GRAT; + case Opm::WellProducer::LRAT: + return LRAT; + case Opm::WellProducer::CRAT: + return CRAT; + case Opm::WellProducer::RESV: + return RESV; + case Opm::WellProducer::BHP: + return BHP; + case Opm::WellProducer::THP: + return THP; + case Opm::WellProducer::GRUP: + return GRUP; + default: + throw std::invalid_argument("unhandled enum value"); + } + } } // namespace ProductionControl @@ -140,6 +156,25 @@ namespace << control << " in input file"); } } + + Mode mode(Opm::WellInjector::ControlModeEnum controlMode) + { + switch ( controlMode ) { + case Opm::WellInjector::GRUP: + return GRUP; + case Opm::WellInjector::RESV: + return RESV; + case Opm::WellInjector::RATE: + return RATE; + case Opm::WellInjector::THP: + return THP; + case Opm::WellInjector::BHP: + return BHP; + default: + throw std::invalid_argument("unhandled enum value"); + } + } + } // namespace InjectionControl @@ -233,6 +268,85 @@ namespace Opm { } + /// Construct wells from deck. + WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState, + const size_t timeStep, + const UnstructuredGrid& grid, + const double* permeability) + : w_(0) + { + if (grid.dimensions != 3) { + OPM_THROW(std::runtime_error, "We cannot initialize wells from a deck unless the corresponding grid is 3-dimensional."); + } + + if (eclipseState->getSchedule()->numWells() == 0) { + OPM_MESSAGE("No wells specified in Schedule section, initializing no wells"); + return; + } + + std::map cartesian_to_compressed; + setupCompressedToCartesian(grid, cartesian_to_compressed); + + // Obtain phase usage data. + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + // These data structures will be filled in this constructor, + // then used to initialize the Wells struct. + std::vector well_names; + std::vector well_data; + + + // For easy lookup: + std::map well_names_to_index; + + ScheduleConstPtr schedule = eclipseState->getSchedule(); + std::vector wells = schedule->getWells(timeStep); + + well_names.reserve(wells.size()); + well_data.reserve(wells.size()); + + createWellsFromSpecs(wells, timeStep, grid, well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability); + setupWellControls(wells, timeStep, well_names, pu); + + { + GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD"); + GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name()); + well_collection_.addField(fieldGroup, timeStep, pu); + addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu); + } + + for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) { + well_collection_.addWell((*wellIter), timeStep, pu); + } + + well_collection_.setWellsPointer(w_); + well_collection_.applyGroupControls(); + + + setupGuideRates(wells, timeStep, well_data, well_names_to_index); + + // Debug output. +#define EXTRA_OUTPUT +#ifdef EXTRA_OUTPUT + /* + std::cout << "\t WELL DATA" << std::endl; + for(int i = 0; i< num_wells; ++i) { + std::cout << i << ": " << well_data[i].type << " " + << well_data[i].control << " " << well_data[i].target + << std::endl; + } + + std::cout << "\n\t PERF DATA" << std::endl; + for(int i=0; i< int(wellperf_data.size()); ++i) { + for(int j=0; j< int(wellperf_data[i].size()); ++j) { + std::cout << i << ": " << wellperf_data[i][j].cell << " " + << wellperf_data[i][j].well_index << std::endl; + } + } + */ +#endif + } + /// Construct wells from deck. @@ -879,4 +993,404 @@ namespace Opm well_collection_.applyExplicitReinjectionControls(well_reservoirrates_phase, well_surfacerates_phase); } + void WellsManager::setupCompressedToCartesian(const UnstructuredGrid& grid, std::map& cartesian_to_compressed ) { + // global_cell is a map from compressed cells to Cartesian grid cells. + // We must make the inverse lookup. + const int* global_cell = grid.global_cell; + + if (global_cell) { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(global_cell[i], i)); + } + } + else { + for (int i = 0; i < grid.number_of_cells; ++i) { + cartesian_to_compressed.insert(std::make_pair(i, i)); + } + } + + } + + void WellsManager::createWellsFromSpecs(std::vector& wells, size_t timeStep, + const UnstructuredGrid& grid, + std::vector& well_names, + std::vector& well_data, + std::map& well_names_to_index, + const PhaseUsage& phaseUsage, + std::map cartesian_to_compressed, + const double* permeability) + { + std::vector > wellperf_data; + wellperf_data.resize(wells.size()); + + int well_index = 0; + for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { + WellConstPtr well = (*wellIter); + { // WELSPECS handling + well_names_to_index[well->name()] = well_index; + well_names.push_back(well->name()); + { + WellData wd; + // If negative (defaulted), set refdepth to a marker + // value, will be changed after getting perforation + // data to the centroid of the cell of the top well + // perforation. + wd.reference_bhp_depth = (well->getRefDepth() < 0.0) ? -1e100 : well->getRefDepth(); + wd.welspecsline = -1; + if (well->isInjector( timeStep )) + wd.type = INJECTOR; + else + wd.type = PRODUCER; + well_data.push_back(wd); + } + } + + { // COMPDAT handling + CompletionSetConstPtr completionSet = well->getCompletions(timeStep); + for (size_t c=0; csize(); c++) { + CompletionConstPtr completion = completionSet->get(c); + int i = completion->getI(); + int j = completion->getJ(); + int k = completion->getK(); + + const int* cpgdim = grid.cartdims; + int cart_grid_indx = i + cpgdim[0]*(j + cpgdim[1]*k); + std::map::const_iterator cgit = cartesian_to_compressed.find(cart_grid_indx); + if (cgit == cartesian_to_compressed.end()) { + OPM_THROW(std::runtime_error, "Cell with i,j,k indices " << i << ' ' << j << ' ' + << k << " not found in grid (well = " << well->name() << ')'); + } + int cell = cgit->second; + PerfData pd; + pd.cell = cell; + if (completion->getCF() > 0.0) { + pd.well_index = completion->getCF(); + } else { + double radius = 0.5*completion->getDiameter(); + if (radius <= 0.0) { + radius = 0.5*unit::feet; + OPM_MESSAGE("**** Warning: Well bore internal radius set to " << radius); + } + std::array cubical = getCubeDim(grid, cell); + const double* cell_perm = &permeability[grid.dimensions*grid.dimensions*cell]; + pd.well_index = computeWellIndex(radius, cubical, cell_perm, completion->getDiameter()); + } + wellperf_data[well_index].push_back(pd); + } + } + well_index++; + } + + // Set up reference depths that were defaulted. Count perfs. + + const int num_wells = well_data.size(); + + int num_perfs = 0; + assert(grid.dimensions == 3); + for (int w = 0; w < num_wells; ++w) { + num_perfs += wellperf_data[w].size(); + if (well_data[w].reference_bhp_depth < 0.0) { + // It was defaulted. Set reference depth to minimum perforation depth. + double min_depth = 1e100; + int num_wperfs = wellperf_data[w].size(); + for (int perf = 0; perf < num_wperfs; ++perf) { + double depth = grid.cell_centroids[3*wellperf_data[w][perf].cell + 2]; + min_depth = std::min(min_depth, depth); + } + well_data[w].reference_bhp_depth = min_depth; + } + } + + // Create the well data structures. + w_ = create_wells(phaseUsage.num_phases, num_wells, num_perfs); + if (!w_) { + OPM_THROW(std::runtime_error, "Failed creating Wells struct."); + } + + + // Add wells. + for (int w = 0; w < num_wells; ++w) { + const int w_num_perf = wellperf_data[w].size(); + std::vector perf_cells(w_num_perf); + std::vector perf_prodind(w_num_perf); + for (int perf = 0; perf < w_num_perf; ++perf) { + perf_cells[perf] = wellperf_data[w][perf].cell; + perf_prodind[perf] = wellperf_data[w][perf].well_index; + } + const double* comp_frac = NULL; + // We initialize all wells with a null component fraction, + // and must (for injection wells) overwrite it later. + int ok = add_well(well_data[w].type, well_data[w].reference_bhp_depth, w_num_perf, + comp_frac, &perf_cells[0], &perf_prodind[0], well_names[w].c_str(), w_); + if (!ok) { + OPM_THROW(std::runtime_error, "Failed adding well " << well_names[w] << " to Wells data structure."); + } + } + + } + + void WellsManager::setupWellControls(std::vector& wells, size_t timeStep, + std::vector& well_names, const PhaseUsage& phaseUsage) { + int well_index = 0; + for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) { + WellConstPtr well = (*wellIter); + + if ( !( well->getStatus( timeStep ) == WellCommon::SHUT || well->getStatus( timeStep ) == WellCommon::OPEN) ) { + OPM_THROW(std::runtime_error, "Currently we do not support well status " << WellCommon::Status2String(well->getStatus( timeStep ))); + } + + if (well->isInjector(timeStep)) { + clear_well_controls(well_index, w_); + int ok = 1; + int control_pos[5] = { -1, -1, -1, -1, -1 }; + + if (well->hasInjectionControl(timeStep , WellInjector::RATE)) { + control_pos[InjectionControl::RATE] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep); + + if (injectorType == WellInjector::TypeEnum::WATER) { + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::OIL) { + distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::GAS) { + distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + + ok = append_well_controls(SURFACE_RATE, + well->getSurfaceInjectionRate( timeStep ) , + distr, + well_index, + w_); + } + + if (ok && well->hasInjectionControl(timeStep , WellInjector::RESV)) { + control_pos[InjectionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + WellInjector::TypeEnum injectorType = well->getInjectorType(timeStep); + + if (injectorType == WellInjector::TypeEnum::WATER) { + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::OIL) { + distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (injectorType == WellInjector::TypeEnum::GAS) { + distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + + ok = append_well_controls(RESERVOIR_RATE, + well->getReservoirInjectionRate( timeStep ), + distr, + well_index, + w_); + } + + if (ok && well->hasInjectionControl(timeStep , WellInjector::BHP)) { + control_pos[InjectionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); + ok = append_well_controls(BHP, + well->getBHPLimit(timeStep), + NULL, + well_index, + w_); + } + + if (ok && well->hasInjectionControl(timeStep , WellInjector::THP)) { + OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]); + } + + + if (!ok) { + OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]); + } + + + { + InjectionControl::Mode mode = InjectionControl::mode( well->getInjectorControlMode(timeStep) ); + int cpos = control_pos[mode]; + if (cpos == -1 && mode != InjectionControl::GRUP) { + OPM_THROW(std::runtime_error, "Control not specified in well " << well_names[well_index]); + } + + // We need to check if the well is shut or not + if (well->getStatus( timeStep ) == WellCommon::SHUT) { + cpos = ~cpos; + } + set_current_control(well_index, cpos, w_); + } + + + // Set well component fraction. + double cf[3] = { 0.0, 0.0, 0.0 }; + if (well->getInjectorType(timeStep) == WellInjector::WATER) { + if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not used, yet found water-injecting well."); + } + cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + } else if (well->getInjectorType(timeStep) == WellInjector::OIL) { + if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-injecting well."); + } + cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + } else if (well->getInjectorType(timeStep) == WellInjector::GAS) { + if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-injecting well."); + } + cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + } + std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases); + + } + + if (well->isProducer(timeStep)) { + // Add all controls that are present in well. + // First we must clear existing controls, in case the + // current WCONPROD line is modifying earlier controls. + clear_well_controls(well_index, w_); + int control_pos[9] = { -1, -1, -1, -1, -1, -1, -1, -1, -1 }; + int ok = 1; + if (ok && well->hasProductionControl(timeStep , WellProducer::ORAT)) { + if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not active and ORAT control specified."); + } + + control_pos[ProductionControl::ORAT] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + ok = append_well_controls(SURFACE_RATE, + -well->getOilRate( timeStep ), + distr, + well_index, + w_); + } + + if (ok && well->hasProductionControl(timeStep , WellProducer::WRAT)) { + if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not active and WRAT control specified."); + } + control_pos[ProductionControl::WRAT] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + ok = append_well_controls(SURFACE_RATE, + -well->getWaterRate(timeStep), + distr, + well_index, + w_); + } + + if (ok && well->hasProductionControl(timeStep , WellProducer::GRAT)) { + if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) { + OPM_THROW(std::runtime_error, "Gas phase not active and GRAT control specified."); + } + control_pos[ProductionControl::GRAT] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0; + ok = append_well_controls(SURFACE_RATE, + -well->getGasRate( timeStep ), + distr, + well_index, + w_); + } + + if (ok && well->hasProductionControl(timeStep , WellProducer::LRAT)) { + if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) { + OPM_THROW(std::runtime_error, "Water phase not active and LRAT control specified."); + } + if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) { + OPM_THROW(std::runtime_error, "Oil phase not active and LRAT control specified."); + } + control_pos[ProductionControl::LRAT] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 0.0, 0.0, 0.0 }; + distr[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0; + distr[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0; + ok = append_well_controls(SURFACE_RATE, + -well->getLiquidRate(timeStep), + distr, + well_index, + w_); + } + + if (ok && well->hasProductionControl(timeStep , WellProducer::RESV)) { + control_pos[ProductionControl::RESV] = well_controls_get_num(w_->ctrls[well_index]); + double distr[3] = { 1.0, 1.0, 1.0 }; + ok = append_well_controls(RESERVOIR_RATE, + -well->getResVRate(timeStep), + distr, + well_index, + w_); + } + + if (ok && well->hasProductionControl(timeStep , WellProducer::BHP)) { + control_pos[ProductionControl::BHP] = well_controls_get_num(w_->ctrls[well_index]); + ok = append_well_controls(BHP, + well->getBHPLimit( timeStep ) , + NULL, + well_index, + w_); + } + + if (ok && well->hasProductionControl(timeStep , WellProducer::THP)) { + OPM_THROW(std::runtime_error, "We cannot handle THP limit for well " << well_names[well_index]); + } + + if (!ok) { + OPM_THROW(std::runtime_error, "Failure occured appending controls for well " << well_names[well_index]); + } + + ProductionControl::Mode mode = ProductionControl::mode(well->getProducerControlMode(timeStep)); + int cpos = control_pos[mode]; + if (cpos == -1 && mode != ProductionControl::GRUP) { + OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]); + } + // If it's shut, we complement the cpos + if (well->getStatus(timeStep) == WellCommon::SHUT) { + cpos = ~cpos; // So we can easily retrieve the cpos later + } + set_current_control(well_index, cpos, w_); + } + well_index++; + } + + } + + void WellsManager::addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage) { + for (auto childIter = parentNode->begin(); childIter != parentNode->end(); ++childIter) { + GroupTreeNodeConstPtr childNode = (*childIter).second; + well_collection_.addGroup(schedule->getGroup(childNode->name()), parentNode->name(), timeStep, phaseUsage); + addChildGroups(childNode, schedule, timeStep, phaseUsage); + } + } + + + + + void WellsManager::setupGuideRates(std::vector& wells, const size_t timeStep, std::vector& well_data, std::map& well_names_to_index) + { + for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) { + WellConstPtr well = *wellIter; + const int wix = well_names_to_index[well->name()]; + WellNode& wellnode = *well_collection_.getLeafNodes()[wix]; + + if (well->getGuideRatePhase(timeStep) != GuideRate::UNDEFINED) { + if (well_data[wix].type == PRODUCER) { + wellnode.prodSpec().guide_rate_ = well->getGuideRate(timeStep); + if (well->getGuideRatePhase(timeStep) == GuideRate::OIL) { + wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL; + } else { + OPM_THROW(std::runtime_error, "Guide rate type " << GuideRate::GuideRatePhaseEnum2String(well->getGuideRatePhase(timeStep)) << " specified for producer " + << well->name() << " in WGRUPCON, cannot handle."); + } + } else if (well_data[wix].type == INJECTOR) { + wellnode.injSpec().guide_rate_ = well->getGuideRate(timeStep); + if (well->getGuideRatePhase(timeStep) == GuideRate::RAT) { + wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT; + } else { + OPM_THROW(std::runtime_error, "Guide rate type " << GuideRate::GuideRatePhaseEnum2String(well->getGuideRatePhase(timeStep)) << " specified for injector " + << well->name() << " in WGRUPCON, cannot handle."); + } + } else { + OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << well->name()); + } + } + } + } + } // namespace Opm diff --git a/opm/core/wells/WellsManager.hpp b/opm/core/wells/WellsManager.hpp index 35f7f471..9e5c314c 100644 --- a/opm/core/wells/WellsManager.hpp +++ b/opm/core/wells/WellsManager.hpp @@ -21,8 +21,11 @@ #define OPM_WELLSMANAGER_HEADER_INCLUDED +#include + #include #include +#include struct Wells; struct UnstructuredGrid; @@ -32,7 +35,22 @@ namespace Opm { class EclipseGridParser; + struct WellData + { + WellType type; + // WellControlType control; + // double target; + double reference_bhp_depth; + // Opm::InjectionSpecification::InjectorType injected_phase; + int welspecsline; + }; + + struct PerfData + { + int cell; + double well_index; + }; /// This class manages a Wells struct in the sense that it /// encapsulates creation and destruction of the wells /// data structure. @@ -40,34 +58,40 @@ namespace Opm class WellsManager { public: - /// Default constructor -- no wells. - WellsManager(); + /// Default constructor -- no wells. + WellsManager(); - /// Construct from existing wells object. - /// WellsManager is not properly initialised in the sense that the logic to - /// manage control switching does not exist. - /// - /// @param[in] W Existing wells object. - WellsManager(struct Wells* W); + /// Construct from existing wells object. + /// WellsManager is not properly initialised in the sense that the logic to + /// manage control switching does not exist. + /// + /// @param[in] W Existing wells object. + WellsManager(struct Wells* W); - /// Construct from input deck and grid. - /// The permeability argument may be zero if the input contain - /// well productivity indices, otherwise it must be given in - /// order to approximate these by the Peaceman formula. - WellsManager(const Opm::EclipseGridParser& deck, - const UnstructuredGrid& grid, - const double* permeability); + /// Construct from input deck and grid. + /// The permeability argument may be zero if the input contain + /// well productivity indices, otherwise it must be given in + /// order to approximate these by the Peaceman formula. + WellsManager(const Opm::EclipseGridParser& deck, + const UnstructuredGrid& grid, + const double* permeability); - /// Destructor. - ~WellsManager(); + + WellsManager(const Opm::EclipseStateConstPtr eclipseState, + const size_t timeStep, + const UnstructuredGrid& grid, + const double* permeability); + + /// Destructor. + ~WellsManager(); /// Does the "deck" define any wells? bool empty() const; - /// Access the managed Wells. - /// The method is named similarly to c_str() in std::string, - /// to make it clear that we are returning a C-compatible struct. - const Wells* c_wells() const; + /// Access the managed Wells. + /// The method is named similarly to c_str() in std::string, + /// to make it clear that we are returning a C-compatible struct. + const Wells* c_wells() const; /// Access the well group hierarchy. const WellCollection& wellCollection() const; @@ -108,18 +132,31 @@ namespace Opm void applyExplicitReinjectionControls(const std::vector& well_reservoirrates_phase, const std::vector& well_surfacerates_phase); + private: - // Disable copying and assignment. - WellsManager(const WellsManager& other); - WellsManager& operator=(const WellsManager& other); + // Disable copying and assignment. + WellsManager(const WellsManager& other); + WellsManager& operator=(const WellsManager& other); + static void setupCompressedToCartesian(const UnstructuredGrid& grid, std::map& cartesian_to_compressed ); + void setupWellControls(std::vector& wells, size_t timeStep, + std::vector& well_names, const PhaseUsage& phaseUsage); - // Data - Wells* w_; + void createWellsFromSpecs( std::vector& wells, size_t timeStep, + const UnstructuredGrid& grid, + std::vector& well_names, + std::vector& well_data, + std::map & well_names_to_index, + const PhaseUsage& phaseUsage, + const std::map cartesian_to_compressed, + const double* permeability); + + void addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage); + void setupGuideRates(std::vector& wells, const size_t timeStep, std::vector& well_data, std::map& well_names_to_index); + + + // Data + Wells* w_; WellCollection well_collection_; - - - - }; } // namespace Opm diff --git a/opm/core/wells/well_controls.c b/opm/core/wells/well_controls.c index a1a0b9b0..5e52d541 100644 --- a/opm/core/wells/well_controls.c +++ b/opm/core/wells/well_controls.c @@ -293,19 +293,38 @@ well_controls_add_new(enum WellControlType type , double target , const double * bool -well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2) +well_controls_equal(const struct WellControls *ctrls1, const struct WellControls *ctrls2 , bool verbose) /* ---------------------------------------------------------------------- */ { bool are_equal = true; - are_equal = (ctrls1->num == ctrls2->num); - are_equal = are_equal && (ctrls1->number_of_phases == ctrls2->number_of_phases); + + if (ctrls1->num != ctrls2->num) { + are_equal = false; + if (verbose) + printf("ctrls1->num:%d ctrls2->num:%d \n",ctrls1->num , ctrls2->num); + } + + if (ctrls1->number_of_phases != ctrls2->number_of_phases) { + are_equal = false; + if (verbose) + printf("ctrls1->number_of_phases:%d ctrls2->number_of_phases:%d \n",ctrls1->number_of_phases , ctrls2->number_of_phases); + } + if (!are_equal) { return are_equal; } - are_equal = are_equal && (memcmp(ctrls1->type, ctrls2->type, ctrls1->num * sizeof *ctrls1->type ) == 0); - are_equal = are_equal && (memcmp(ctrls1->target, ctrls2->target, ctrls1->num * sizeof *ctrls1->target ) == 0); - are_equal = are_equal && (memcmp(ctrls1->distr, ctrls2->distr, ctrls1->num * ctrls1->number_of_phases * sizeof *ctrls1->distr ) == 0); + if (memcmp(ctrls1->type, ctrls2->type, ctrls1->num * sizeof *ctrls1->type ) != 0) { + are_equal = false; + if (verbose) + printf("The ->type vectors are different \n"); + } + + if (memcmp(ctrls1->target, ctrls2->target, ctrls1->num * sizeof *ctrls1->target ) != 0) { + are_equal = false; + if (verbose) + printf("The ->target vectors are different \n"); + } are_equal = are_equal && (ctrls1->cpty == ctrls2->cpty); return are_equal; diff --git a/opm/core/wells/wells.c b/opm/core/wells/wells.c index 13b8c69b..81a8ddc1 100644 --- a/opm/core/wells/wells.c +++ b/opm/core/wells/wells.c @@ -541,7 +541,7 @@ clone_wells(const struct Wells *W) /* ---------------------------------------------------------------------- */ bool -wells_equal(const struct Wells *W1, const struct Wells *W2) +wells_equal(const struct Wells *W1, const struct Wells *W2 , bool verbose) /* ---------------------------------------------------------------------- */ { bool are_equal = true; @@ -567,10 +567,25 @@ wells_equal(const struct Wells *W1, const struct Wells *W2) are_equal = are_equal && (strcmp(W1->name[i], W2->name[i]) == 0); else are_equal = are_equal && (W1->name[i] == W2->name[i]); + + if (verbose && !are_equal) + printf("Well name[%d] %s and %s are different \n", i , W1->name[i] , W2->name[i]); + } + if (W1->type[i] != W2->type[i]) { + are_equal = false; + if (verbose) + printf("Well->type[%d] different %d %d \n",i , W1->type[i] , W2->type[i] ); + } + if (W1->depth_ref[i] != W2->depth_ref[i]) { + are_equal = false; + if (verbose) + printf("Well->depth_ref[%d] different %g %g \n",i , W1->depth_ref[i] , W2->depth_ref[i] ); + } + if (!well_controls_equal(W1->ctrls[i], W2->ctrls[i],verbose)) { + are_equal = false; + if (verbose) + printf("Well controls are different for well[%d]:%s \n",i,W1->name[i]); } - are_equal = are_equal && (W1->type[i] == W2->type[i]); - are_equal = are_equal && (W1->depth_ref[i] == W2->depth_ref[i]); - are_equal = are_equal && (well_controls_equal(W1->ctrls[i], W2->ctrls[i])); } @@ -581,10 +596,16 @@ wells_equal(const struct Wells *W1, const struct Wells *W2) are_equal = are_equal && (mgmt1->well_cpty == mgmt2->well_cpty); } - are_equal = are_equal && (memcmp(W1->comp_frac, W2->comp_frac, W1->number_of_wells * W1->number_of_phases * sizeof *W1->comp_frac ) == 0); - are_equal = are_equal && (memcmp(W1->well_connpos, W2->well_connpos, (1 + W1->number_of_wells) * sizeof *W1->well_connpos ) == 0); - if (!are_equal) { - return are_equal; + if (memcmp(W1->comp_frac, W2->comp_frac, W1->number_of_wells * W1->number_of_phases * sizeof *W1->comp_frac ) != 0) { + are_equal = false; + if (verbose) + printf("Component fractions different \n"); + } + + if (memcmp(W1->well_connpos, W2->well_connpos, (1 + W1->number_of_wells) * sizeof *W1->well_connpos ) != 0) { + are_equal = false; + if (verbose) + printf("perforation position map difference \n"); } { diff --git a/tests/testBlackoilState1.DATA b/tests/testBlackoilState1.DATA index 5bbb4eb0..4220b5b4 100644 --- a/tests/testBlackoilState1.DATA +++ b/tests/testBlackoilState1.DATA @@ -10,5 +10,9 @@ DYV DZV 10*10 / +ACTNUM + 1 998*2 3 / + + DEPTHZ 121*2000 / diff --git a/tests/test_parser.cpp b/tests/test_parser.cpp new file mode 100644 index 00000000..ed0c1cd3 --- /dev/null +++ b/tests/test_parser.cpp @@ -0,0 +1,54 @@ +/* + 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 + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE OPM-ParserTest +#include + +#include +#include +#include + +#include +#include +#include +#include + +BOOST_AUTO_TEST_CASE(CreateParser) +{ + const std::string filename1 = "testBlackoilState1.DATA"; + Opm::ParserPtr parser(new Opm::Parser() ); + Opm::DeckConstPtr deck = parser->parseFile( filename1 ); + + BOOST_CHECK_EQUAL( 6U , deck->size() ); + Opm::DeckItemConstPtr actnum = deck->getKeyword("ACTNUM")->getRecord(0)->getItem(0); + const std::vector& actnum_data = actnum->getIntData(); + + BOOST_CHECK_EQUAL( 1000U , actnum->size() ); + BOOST_CHECK_EQUAL( 1, actnum_data[0] ); + BOOST_CHECK_EQUAL( 2, actnum_data[400] ); + BOOST_CHECK_EQUAL( 3, actnum_data[999] ); +} diff --git a/tests/test_wellcollection.cpp b/tests/test_wellcollection.cpp new file mode 100644 index 00000000..d22ba3b4 --- /dev/null +++ b/tests/test_wellcollection.cpp @@ -0,0 +1,91 @@ +/* + 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 + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellCollectionTest +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace Opm; + +BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) { + ParserPtr parser(new Parser()); + boost::filesystem::path scheduleFile("wells_group.data"); + DeckConstPtr deck = parser->parseFile(scheduleFile.string()); + EclipseStateConstPtr eclipseState(new EclipseState(deck)); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + GroupTreeNodePtr field=eclipseState->getSchedule()->getGroupTree(2)->getNode("FIELD"); + GroupTreeNodePtr g1=eclipseState->getSchedule()->getGroupTree(2)->getNode("G1"); + GroupTreeNodePtr g2=eclipseState->getSchedule()->getGroupTree(2)->getNode("G2"); + + WellCollection collection; + + // Add groups to WellCollection + GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(field->name()); + collection.addField(fieldGroup, 2, pu); + + for (auto iter = field->begin(); iter != field->end(); ++iter) { + GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name()); + collection.addGroup(childGroupNode, fieldGroup->name(), 2, pu); + } + + GroupConstPtr g1Group = eclipseState->getSchedule()->getGroup(g1->name()); + for (auto iter = g1->begin(); iter != g1->end(); ++iter) { + GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name()); + collection.addGroup(childGroupNode, g1Group->name(), 2, pu); + } + + + GroupConstPtr g2Group = eclipseState->getSchedule()->getGroup(g2->name()); + for (auto iter = g2->begin(); iter != g2->end(); ++iter) { + GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name()); + collection.addGroup(childGroupNode, g2Group->name(), 2, pu); + } + + BOOST_CHECK_EQUAL("FIELD", collection.findNode("FIELD")->name()); + BOOST_CHECK_EQUAL("FIELD", collection.findNode("G1")->getParent()->name()); + BOOST_CHECK_EQUAL("FIELD", collection.findNode("G2")->getParent()->name()); + + // Add wells to WellCollection + WellCollection wellCollection; + std::vector wells = eclipseState->getSchedule()->getWells(); + for (size_t i=0; igetParent()->name()); + BOOST_CHECK_EQUAL("G1", collection.findNode("INJ2")->getParent()->name()); + BOOST_CHECK_EQUAL("G2", collection.findNode("PROD1")->getParent()->name()); + BOOST_CHECK_EQUAL("G2", collection.findNode("PROD2")->getParent()->name()); +} + diff --git a/tests/test_wells.cpp b/tests/test_wells.cpp index 9e76d899..8502c4b8 100644 --- a/tests/test_wells.cpp +++ b/tests/test_wells.cpp @@ -204,7 +204,7 @@ BOOST_AUTO_TEST_CASE(Copy) BOOST_CHECK_EQUAL( dist1[p] , dist2[p]); } } - BOOST_CHECK( well_controls_equal( c1 , c2 )); + BOOST_CHECK( well_controls_equal( c1 , c2 , false) ); } } } @@ -219,7 +219,7 @@ BOOST_AUTO_TEST_CASE(Equals_WellsEqual_ReturnsTrue) { std::shared_ptr W2(create_wells(nphases, nwells, nperfs), destroy_wells); - BOOST_CHECK(wells_equal(W1.get(), W2.get())); + BOOST_CHECK(wells_equal(W1.get(), W2.get() , false)); } @@ -232,5 +232,5 @@ BOOST_AUTO_TEST_CASE(Equals_WellsDiffer_ReturnsFalse) { std::shared_ptr W2(create_wells(nphases, 3, nperfs), destroy_wells); - BOOST_CHECK(!wells_equal(W1.get(), W2.get())); + BOOST_CHECK(!wells_equal(W1.get(), W2.get() , false )); } diff --git a/tests/test_wellsgroup.cpp b/tests/test_wellsgroup.cpp new file mode 100644 index 00000000..9ba9e0ef --- /dev/null +++ b/tests/test_wellsgroup.cpp @@ -0,0 +1,108 @@ +/* + Copyright 2014 Statoil. + + 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 + +#if HAVE_DYNAMIC_BOOST_TEST +#define BOOST_TEST_DYN_LINK +#endif + +#define NVERBOSE // Suppress own messages when throw()ing + +#define BOOST_TEST_MODULE WellsGroupTest + +#include +#include + +#include +#include + +#include +#include + +#include +#include +#include + +#include +#include + +using namespace Opm; + +BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) { + ParserPtr parser(new Parser()); + boost::filesystem::path scheduleFile("wells_group.data"); + DeckConstPtr deck = parser->parseFile(scheduleFile.string()); + EclipseStateConstPtr eclipseState(new EclipseState(deck)); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + std::vector wells = eclipseState->getSchedule()->getWells(); + + for (size_t i=0; i wellsGroup = createWellWellsGroup(well, 2, pu); + BOOST_CHECK_EQUAL(well->name(), wellsGroup->name()); + if (well->isInjector(2)) { + BOOST_CHECK_EQUAL(well->getSurfaceInjectionRate(2), wellsGroup->injSpec().surface_flow_max_rate_); + BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->injSpec().BHP_limit_); + BOOST_CHECK_EQUAL(well->getReservoirInjectionRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_); + } + if (well->isProducer(2)) { + BOOST_CHECK_EQUAL(well->getResVRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->prodSpec().BHP_limit_); + BOOST_CHECK_EQUAL(well->getOilRate(2), wellsGroup->prodSpec().oil_max_rate_); + BOOST_CHECK_EQUAL(well->getWaterRate(2), wellsGroup->prodSpec().water_max_rate_); + BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_); + } + } +} + + +BOOST_AUTO_TEST_CASE(ConstructGroupFromGroup) { + ParserPtr parser(new Parser()); + boost::filesystem::path scheduleFile("wells_group.data"); + DeckConstPtr deck = parser->parseFile(scheduleFile.string()); + EclipseStateConstPtr eclipseState(new EclipseState(deck)); + PhaseUsage pu = phaseUsageFromDeck(eclipseState); + + std::vector nodes = eclipseState->getSchedule()->getGroupTree(2)->getNodes(); + + for (size_t i=0; igetSchedule()->getGroup(nodes[i]->name()); + std::shared_ptr wellsGroup = createGroupWellsGroup(group, 2, pu); + BOOST_CHECK_EQUAL(group->name(), wellsGroup->name()); + if (group->isInjectionGroup(2)) { + BOOST_CHECK_EQUAL(group->getSurfaceMaxRate(2), wellsGroup->injSpec().surface_flow_max_rate_); + BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(group->getTargetReinjectFraction(2), wellsGroup->injSpec().reinjection_fraction_target_); + BOOST_CHECK_EQUAL(group->getTargetVoidReplacementFraction(2), wellsGroup->injSpec().voidage_replacment_fraction_); + } + if (group->isProductionGroup(2)) { + BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_); + BOOST_CHECK_EQUAL(group->getGasTargetRate(2), wellsGroup->prodSpec().gas_max_rate_); + BOOST_CHECK_EQUAL(group->getOilTargetRate(2), wellsGroup->prodSpec().oil_max_rate_); + BOOST_CHECK_EQUAL(group->getWaterTargetRate(2), wellsGroup->prodSpec().water_max_rate_); + BOOST_CHECK_EQUAL(group->getLiquidTargetRate(2), wellsGroup->prodSpec().liquid_max_rate_); + } + } +} + + + diff --git a/tests/test_wellsmanager.cpp b/tests/test_wellsmanager.cpp index b3465fac..8cb097e2 100644 --- a/tests/test_wellsmanager.cpp +++ b/tests/test_wellsmanager.cpp @@ -27,6 +27,10 @@ #define BOOST_TEST_MODULE WellsManagerTests #include + +#include +#include + #include #include #include @@ -57,7 +61,7 @@ void wells_static_check(const Wells* wells) { } -/* +/* The number of controls is determined by looking at which elements have been given explicit - non-default - values in the WCONxxxx keyword. Is that at all interesting? @@ -67,7 +71,7 @@ void wells_static_check(const Wells* wells) { void check_controls_epoch0( struct WellControls ** ctrls) { // The injector { - const struct WellControls * ctrls0 = ctrls[0]; + const struct WellControls * ctrls0 = ctrls[0]; BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0)); // The number of controls for the injector == 3?? BOOST_CHECK_EQUAL( SURFACE_RATE , well_controls_iget_type(ctrls0 , 0) ); @@ -83,7 +87,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) { BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls0) ); // The phase distribution in the active target - { + { const double * distr = well_controls_iget_distr( ctrls0 , 0 ); BOOST_CHECK_EQUAL( 0 , distr[0] ); // Water BOOST_CHECK_EQUAL( 0 , distr[1] ); // Oil @@ -106,7 +110,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) { BOOST_CHECK_EQUAL( 0 , well_controls_get_current(ctrls1)); // The phase distribution in the active target - { + { const double * distr = well_controls_iget_distr( ctrls1 , 0 ); BOOST_CHECK_EQUAL( 0 , distr[0] ); // Water BOOST_CHECK_EQUAL( 1 , distr[1] ); // Oil @@ -121,7 +125,7 @@ void check_controls_epoch0( struct WellControls ** ctrls) { void check_controls_epoch1( struct WellControls ** ctrls) { // The injector { - const struct WellControls * ctrls0 = ctrls[0]; + const struct WellControls * ctrls0 = ctrls[0]; BOOST_CHECK_EQUAL( 3 , well_controls_get_num(ctrls0)); // The number of controls for the injector == 3?? BOOST_CHECK_EQUAL( SURFACE_RATE , well_controls_iget_type(ctrls0 , 0 )); @@ -136,7 +140,7 @@ void check_controls_epoch1( struct WellControls ** ctrls) { // Which control is active BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls0)); - { + { const double * distr = well_controls_iget_distr( ctrls0 , 1 ); BOOST_CHECK_EQUAL( 1 , distr[0] ); // Water BOOST_CHECK_EQUAL( 0 , distr[1] ); // Oil @@ -160,7 +164,7 @@ void check_controls_epoch1( struct WellControls ** ctrls) { // Which control is active BOOST_CHECK_EQUAL( 1 , well_controls_get_current(ctrls1) ); - { + { const double * distr = well_controls_iget_distr( ctrls1 , 1 ); BOOST_CHECK_EQUAL( 1 , distr[0] ); // Water BOOST_CHECK_EQUAL( 1 , distr[1] ); // Oil @@ -170,8 +174,6 @@ void check_controls_epoch1( struct WellControls ** ctrls) { } - - BOOST_AUTO_TEST_CASE(Constructor_Works) { Opm::EclipseGridParser Deck("wells_manager_data.data"); Opm::GridManager gridManager(Deck); @@ -190,11 +192,86 @@ BOOST_AUTO_TEST_CASE(Constructor_Works) { Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL); const Wells* wells = wellsManager.c_wells(); - wells_static_check( wells ); + wells_static_check( wells ); check_controls_epoch1( wells->ctrls ); } } +BOOST_AUTO_TEST_CASE(New_Constructor_Works) { + + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data"))); + + Opm::EclipseGridParser Deck("wells_manager_data.data"); + Opm::GridManager gridManager(Deck); + + Deck.setCurrentEpoch(0); + { + Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL); + Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL); + + std::cout << "Checking new well structure, epoch 0" << std::endl; + wells_static_check( wellsManager.c_wells() ); + + std::cout << "Checking old well structure, epoch 0" << std::endl; + wells_static_check( oldWellsManager.c_wells() ); + + check_controls_epoch0( wellsManager.c_wells()->ctrls ); + + BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells() , false)); + } + + Deck.setCurrentEpoch(1); + { + Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL); + Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL); + + std::cout << "Checking new well structure, epoch 1" << std::endl; + wells_static_check( wellsManager.c_wells() ); + + std::cout << "Checking old well structure, epoch 1" << std::endl; + wells_static_check( oldWellsManager.c_wells() ); + + check_controls_epoch1( wellsManager.c_wells()->ctrls ); + + BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(),false)); + } +} + + +BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) { + + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_expanded.data"))); + + Opm::EclipseGridParser Deck("wells_manager_data_expanded.data"); + Opm::GridManager gridManager(Deck); + + Deck.setCurrentEpoch(0); + { + Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL); + Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL); + + BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells(),false)); + } + + Deck.setCurrentEpoch(1); + { + Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL); + Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL); + + BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true)); + } + + Deck.setCurrentEpoch(2); + { + Opm::WellsManager wellsManager(eclipseState, 2, *gridManager.c_grid(), NULL); + Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL); + + BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true)); + } +} + BOOST_AUTO_TEST_CASE(WellsEqual) { @@ -207,8 +284,8 @@ BOOST_AUTO_TEST_CASE(WellsEqual) { Deck.setCurrentEpoch(1); Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL); - BOOST_CHECK( wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells()) ); - BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells()) ); + BOOST_CHECK( wells_equal( wellsManager0.c_wells() , wellsManager0.c_wells(),false) ); + BOOST_CHECK( !wells_equal( wellsManager0.c_wells() , wellsManager1.c_wells(),false) ); } @@ -222,15 +299,30 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) { Deck.setCurrentEpoch(1); Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL); - BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0])); - BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1])); - BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0])); - BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1])); + BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false)); + BOOST_CHECK( well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false)); + BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager1.c_wells()->ctrls[0] , false)); + BOOST_CHECK( well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager1.c_wells()->ctrls[1] , false)); - BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1])); - BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0])); - BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0])); - BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1])); + BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[1] , false)); + BOOST_CHECK( !well_controls_equal( wellsManager0.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[0] , false)); + BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[0] , wellsManager0.c_wells()->ctrls[0] , false)); + BOOST_CHECK( !well_controls_equal( wellsManager1.c_wells()->ctrls[1] , wellsManager0.c_wells()->ctrls[1] , false)); } + +BOOST_AUTO_TEST_CASE(WellHasSTOP_ExceptionIsThrown) { + Opm::EclipseGridParser Deck("wells_manager_data_wellSTOP.data"); + Opm::GridManager gridManager(Deck); + + Opm::ParserPtr parser(new Opm::Parser()); + Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_wellSTOP.data"))); + + Deck.setCurrentEpoch(0); + + BOOST_CHECK_THROW( new Opm::WellsManager(eclipseState, 0, *gridManager.c_grid(), NULL), std::runtime_error ); +} + + + diff --git a/tests/wells_group.data b/tests/wells_group.data new file mode 100755 index 00000000..5d2e58cb --- /dev/null +++ b/tests/wells_group.data @@ -0,0 +1,68 @@ +OIL +GAS +WATER + +DIMENS + 10 10 5 / + +GRID + +DXV +10*1000.0 / + +DYV +10*1000.0 / + +DZV +10.0 20.0 30.0 10.0 5.0 / + +DEPTHZ +121*2000 +/ + +SCHEDULE + +GRUPTREE + 'G1' 'FIELD' / + 'G2' 'FIELD' / +/ + + +WELSPECS + 'INJ1' 'G1' 1 1 8335 'GAS' / + 'PROD1' 'G2' 10 10 8400 'OIL' / +/ + +TSTEP + 14.0 / +/ + +WELSPECS + 'INJ2' 'G1' 1 1 8335 'GAS' / + 'PROD2' 'G2' 10 10 8400 'OIL' / +/ + +GCONINJE +'G1' GAS RATE 30000 / +/ + +GCONPROD +'G2' ORAT 10000 / +/ + +WCONINJE + 'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 / + 'INJ2' 'WATER' 'OPEN' 'RESV' 10 20 40 / +/ + +WCONPROD + 'PROD1' 'OPEN' 'RESV' 999 3* 123 100 / + 'PROD2' 'OPEN' 'RESV' 999 3* 123 100 / +/ + +TSTEP + 3 / +/ + + +END diff --git a/tests/wells_manager_data.data b/tests/wells_manager_data.data index f22152c1..732b528b 100755 --- a/tests/wells_manager_data.data +++ b/tests/wells_manager_data.data @@ -17,8 +17,8 @@ DZV 10.0 20.0 30.0 10.0 5.0 / DEPTHZ -121*2000 -/ +121*2000 / + SCHEDULE @@ -34,11 +34,11 @@ COMPDAT WCONPROD 'PROD1' 'OPEN' 'ORAT' 20000 4* 1000 / - / +/ WCONINJE 'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 / - / +/ DATES @@ -53,10 +53,26 @@ WCONINJE 'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 / / - +END TSTEP -14.0 - / + 14.0 / +/ + + +WELSPECS + 'TEST1' 'G1' 1 1 8335 'GAS' / + 'TEST2' 'G2' 10 10 8400 'OIL' / +/ + +GRUPTREE + 'G1' 'SHIP' / + 'G2' 'SHIP' / +/ + +TSTEP + 3 / +/ + END diff --git a/tests/wells_manager_data_expanded.data b/tests/wells_manager_data_expanded.data new file mode 100755 index 00000000..67207768 --- /dev/null +++ b/tests/wells_manager_data_expanded.data @@ -0,0 +1,86 @@ +OIL +GAS +WATER + +DIMENS + 10 10 5 / + +GRID + +DXV +10*1000.0 / + +DYV +10*1000.0 / + +DZV +10.0 20.0 30.0 10.0 5.0 / + +-- The DEPTHZ keyword is only here to satisfy the old parser; content might +-- completely bogus. +DEPTHZ +121*2000 / + + +SCHEDULE + +WELSPECS + 'INJ1' 'G' 1 1 8335 'GAS' / + 'PROD1' 'G' 10 10 8400 'OIL' / +/ + +COMPDAT + 'INJ1' 1 1 1 2 'OPEN' 1 10.6092 0.5 / + 'INJ1' 1 1 3 5 'OPEN' 1 12.6092 0.5 / + 'INJ1' 2 2 1 1 'OPEN' 1 14.6092 0.5 / + 'PROD1' 10 3 3 3 'OPEN' 0 10.6092 0.5 / +/ + +WCONPROD + 'PROD1' 'OPEN' 'ORAT' 20000 4* 1000 / +/ + +WCONINJE + 'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 / +/ + + +DATES + 1 'FEB' 2000 / +/ + +WELSPECS + 'INJ1' 'G3' 1 1 8335 'GAS' / + 'QNJ2' 'G3' 1 1 8335 'GAS' / +/ + + +COMPDAT + 'QNJ2' 3 4 1 2 'OPEN' 1 10.6092 0.5 / + 'QNJ2' 4 4 3 5 'OPEN' 1 12.6092 0.5 / +/ + +WCONPROD + 'PROD1' 'OPEN' 'RESV' 999 3* 123 100 / +/ + +WCONINJE + 'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 / + 'QNJ2' 'WATER' 'OPEN' 'RESV' 7 33 39 / +/ + + + +TSTEP +14.0 / +/ + +WELOPEN + 'INJ1' 'SHUT' 5* / +/ + +TSTEP +14.0 / +/ + +END diff --git a/tests/wells_manager_data_wellSTOP.data b/tests/wells_manager_data_wellSTOP.data new file mode 100755 index 00000000..bdcbc381 --- /dev/null +++ b/tests/wells_manager_data_wellSTOP.data @@ -0,0 +1,42 @@ +OIL +GAS +WATER + +DIMENS + 10 10 5 / + +GRID + +DXV +10*1000.0 / + +DYV +10*1000.0 / + +DZV +10.0 20.0 30.0 10.0 5.0 / + +DEPTHZ +121*2000 +/ + +SCHEDULE + +WELSPECS + 'INJ1' 'G' 1 1 8335 'GAS' / + 'PROD1' 'G' 10 10 8400 'OIL' / +/ + +COMPDAT + 'INJ1' 1 1 1 1 'OPEN' 1 10.6092 0.5 / + 'PROD1' 10 3 3 3 'OPEN' 0 10.6092 0.5 / +/ + +WELOPEN + 'INJ1' 'STOP' 5* / +/ + +TSTEP + 10 / + +END