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

Merge branch 'master' of git://github.com/OPM/opm-core

Browse Source 
Conflicts:
	CMakeLists_files.cmake
This commit is contained in:
osae 2014-03-26 14:51:06 +01:00
commit 5df819d509
88 changed files with 4443 additions and 1060 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
CMakeLists.txt
View File

@ -46,6 +46,7 @@ include (CMakeLists_files.cmake)
macro (config_hook)
# opm_need_version_of ("dune-common")
opm_need_version_of ("dune-istl")
list (APPEND ${project}_CONFIG_IMPL_VARS
HAVE_DUNE_ISTL
)

23
CMakeLists_files.cmake
View File

@ -162,15 +162,19 @@ list (APPEND TEST_SOURCE_FILES
tests/test_wachspresscoord.cpp
tests/test_column_extract.cpp
tests/test_geom2d.cpp
tests/test_linearsolver.cpp
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_equil.cpp
tests/test_regionmapping.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:
@ -179,14 +183,17 @@ 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/capillary.DATA
tests/capillary_overlap.DATA
tests/deadfluids.DATA
tests/equil_liveoil.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:

7
cmake/Modules/Findcjson.cmake
View File

@ -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"

36
cmake/Modules/Findopm-parser.cmake
View File

@ -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 ()

1
cmake/Modules/OpmFind.cmake
View File

@ -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

5
cmake/Modules/UseDynamicBoost.cmake
View File

@ -1,5 +1,6 @@
# be sure that this component is searched for
find_package (Boost COMPONENTS unit_test_framework QUIET)
if (NOT ${Boost_UNIT_TEST_FRAMEWORK_FOUND})
find_package (Boost 1.44.0 COMPONENTS unit_test_framework QUIET)
endif ()
if (${Boost_UNIT_TEST_FRAMEWORK_FOUND})
# setup to do a test compile

4
cmake/Modules/dune-cornerpoint-prereqs.cmake
View File

@ -15,8 +15,8 @@ set (dune-cornerpoint_DEPS
# compile with C++0x/11 support if available
"CXX11Features"
# 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 iostreams unit_test_framework REQUIRED"
# DUNE dependency
"dune-common REQUIRED;
dune-grid REQUIRED;

4
cmake/Modules/opm-autodiff-prereqs.cmake
View File

@ -12,8 +12,8 @@ set (opm-autodiff_DEPS
# Compile with C++0x/11 support if available
"CXX11Features"
# 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 iostreams unit_test_framework REQUIRED"
# DUNE prerequisites
"dune-common REQUIRED;
dune-istl REQUIRED;

18
cmake/Modules/opm-benchmarks-prereqs.cmake Normal file
View File

@ -0,0 +1,18 @@
# -*- mode: cmake; tab-width: 2; indent-tabs-mode: t; truncate-lines: t; compile-command: "cmake -Wdev" -*-
# vim: set filetype=cmake autoindent tabstop=2 shiftwidth=2 noexpandtab softtabstop=2 nowrap:
# defines that must be present in config.h for our headers
set (opm-benchmarks_CONFIG_VAR
)
# dependencies
set (opm-benchmarks_DEPS
# compile with C++0x/11 support if available
"CXX11Features REQUIRED"
# various runtime library enhancements
"Boost 1.44.0
COMPONENTS date_time filesystem system iostreams unit_test_framework REQUIRED"
# OPM dependency
"opm-core"
"opm-upscaling"
)

6
cmake/Modules/opm-core-prereqs.cmake
View File

@ -14,8 +14,8 @@ set (opm-core_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 iostreams unit_test_framework REQUIRED"
# matrix library
"BLAS REQUIRED"
"LAPACK REQUIRED"
@ -30,4 +30,6 @@ set (opm-core_DEPS
# DUNE dependency
"dune-common"
"dune-istl"
# Parser library for ECL-type simulation models
"opm-parser REQUIRED"
)

21
cmake/Modules/opm-parser-prereqs.cmake
View File

@ -3,16 +3,17 @@
# defines that must be present in config.h for our headers
set (opm-parser_CONFIG_VAR
HAVE_ERT
)
HAVE_ERT
)
# dependencies
set (opm-parser_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"
)
# 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 iostream unit_test_framework REQUIRED"
"cJSON"
)

4
cmake/Modules/opm-polymer-prereqs.cmake
View File

@ -12,8 +12,8 @@ set (opm-polymer_DEPS
# compile with C++0x/11 support if available
"CXX11Features"
# 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 iostreams unit_test_framework REQUIRED"
# Ensembles-based Reservoir Tools
"ERT"
# OPM dependency

4
cmake/Modules/opm-porsol-prereqs.cmake
View File

@ -12,8 +12,8 @@ set (opm-porsol_DEPS
# compile with C++0x/11 support if available
"CXX11Features"
# 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 iostreams unit_test_framework REQUIRED"
# DUNE dependency
"dune-common REQUIRED;
dune-istl REQUIRED;

4
cmake/Modules/opm-upscaling-prereqs.cmake
View File

@ -13,8 +13,8 @@ set (opm-upscaling_DEPS
# compile with C++0x/11 support if available
"CXX11Features"
# 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 iostreams unit_test_framework REQUIRED"
# matrix library
"BLAS REQUIRED"
"LAPACK REQUIRED"

4
cmake/Modules/opm-verteq-prereqs.cmake
View File

@ -12,8 +12,8 @@ set (opm-verteq_DEPS
# compile with C++0x/11 support if available
"CXX11Features"
# 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 iostreams unit_test_framework REQUIRED"
# OPM dependency
"opm-core"
)

2
dune.module
View File

@ -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

8
examples/compute_tof.cpp
View File

@ -47,6 +47,9 @@
#include <opm/core/tof/TofReorder.hpp>
#include <opm/core/tof/TofDiscGalReorder.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <memory>
#include <boost/filesystem.hpp>
@ -111,13 +114,16 @@ try
double gravity[3] = { 0.0 };
if (use_deck) {
std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile(deck_filename)));
deck.reset(new EclipseGridParser(deck_filename));
// Grid init
grid.reset(new GridManager(*deck));
// Rock and fluid init
props.reset(new IncompPropertiesFromDeck(*deck, *grid->c_grid()));
// Wells init.
wells.reset(new Opm::WellsManager(*deck, *grid->c_grid(), props->permeability()));
wells.reset(new Opm::WellsManager(eclipseState , 0 , *grid->c_grid(), props->permeability()));
// Gravity.
gravity[2] = deck->hasField("NOGRAV") ? 0.0 : unit::gravity;
// Init state variables (saturation and pressure).

10
examples/sim_2p_comp_reorder.cpp
View File

@ -44,6 +44,10 @@
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/SimulatorCompressibleTwophase.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <memory>
#include <boost/filesystem.hpp>
@ -80,6 +84,7 @@ try
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
EclipseStateConstPtr eclipseState;
std::unique_ptr<EclipseGridParser> deck;
std::unique_ptr<GridManager> grid;
std::unique_ptr<BlackoilPropertiesInterface> props;
@ -89,7 +94,10 @@ try
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
ParserPtr parser(new Opm::Parser());
std::string deck_filename = param.get<std::string>("deck_filename");
eclipseState.reset( new EclipseState(parser->parseFile(deck_filename)));
deck.reset(new EclipseGridParser(deck_filename));
// Grid init
grid.reset(new GridManager(*deck));
@ -242,7 +250,7 @@ try
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.

10
examples/sim_2p_incomp.cpp
View File

@ -45,6 +45,10 @@
#include <opm/core/simulator/WellState.hpp>
#include <opm/core/simulator/SimulatorIncompTwophase.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <memory>
#include <boost/filesystem.hpp>
@ -94,6 +98,8 @@ try
// If we have a "deck_filename", grid and props will be read from that.
bool use_deck = param.has("deck_filename");
EclipseStateConstPtr eclipseState;
std::unique_ptr<EclipseGridParser> deck;
std::unique_ptr<GridManager> grid;
std::unique_ptr<IncompPropertiesInterface> props;
@ -103,8 +109,10 @@ try
// int max_well_control_iterations = 0;
double gravity[3] = { 0.0 };
if (use_deck) {
ParserPtr parser(new Opm::Parser());
std::string deck_filename = param.get<std::string>("deck_filename");
deck.reset(new EclipseGridParser(deck_filename));
eclipseState.reset( new EclipseState(parser->parseFile(deck_filename)));
// Grid init
grid.reset(new GridManager(*deck));
// Rock and fluid init
@ -260,7 +268,7 @@ try
<< simtimer.numSteps() - step << ")\n\n" << std::flush;
// Create new wells, well_state
WellsManager wells(*deck, *grid->c_grid(), props->permeability());
WellsManager wells(eclipseState , epoch , *grid->c_grid(), props->permeability());
// @@@ HACK: we should really make a new well state and
// properly transfer old well state to it every epoch,
// since number of wells may change etc.

7
examples/wells_example.cpp
View File

@ -19,6 +19,9 @@
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/props/rock/RockCompressibility.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
int main(int argc, char** argv)
try
{
@ -39,9 +42,11 @@ try
// Define rock and fluid properties
IncompPropertiesFromDeck incomp_properties(parser, *grid.c_grid());
RockCompressibility rock_comp(parser);
ParserPtr newParser(new Opm::Parser());
EclipseStateConstPtr eclipseState(new Opm::EclipseState(newParser->parseFile(file_name)));
// Finally handle the wells
WellsManager wells(parser, *grid.c_grid(), incomp_properties.permeability());
WellsManager wells(eclipseState , 0 , *grid.c_grid(), incomp_properties.permeability());
double gravity[3] = {0.0, 0.0, parameters.getDefault<double>("gravity", 0.0)};
Opm::LinearSolverFactory linsolver(parameters);

161
opm/core/grid/GridManager.cpp
View File

@ -23,16 +23,13 @@
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/cornerpoint_grid.h>
#include <array>
#include <algorithm>
#include <numeric>
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<double>& zcorn = newParserDeck->getKeyword("ZCORN")->getSIDoubleData();
const std::vector<double>& coord = newParserDeck->getKeyword("COORD")->getSIDoubleData();
const int* actnum = NULL;
if (newParserDeck->hasKeyword("ACTNUM")) {
actnum = &(newParserDeck->getKeyword("ACTNUM")->getIntData()[0]);
}
std::array<int, 3> 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<double*>(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<int, 3> 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<double>& dxv = newParserDeck->getKeyword("DXV")->getSIDoubleData();
const std::vector<double>& dyv = newParserDeck->getKeyword("DYV")->getSIDoubleData();
const std::vector<double>& dzv = newParserDeck->getKeyword("DZV")->getSIDoubleData();
std::vector<double> x = coordsFromDeltas(dxv);
std::vector<double> y = coordsFromDeltas(dyv);
std::vector<double> 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<double> top_depths_vec;
if (newParserDeck->hasKeyword("DEPTHZ")) {
const std::vector<double>& 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<double>& 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

15
opm/core/grid/GridManager.hpp
View File

@ -20,10 +20,12 @@
#ifndef OPM_GRIDMANAGER_HEADER_INCLUDED
#define OPM_GRIDMANAGER_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <string>
struct UnstructuredGrid;
struct grdecl;
namespace Opm
{
@ -42,7 +44,10 @@ namespace Opm
{
public:
/// Construct a 3d corner-point grid or tensor grid from a deck.
GridManager(const Opm::EclipseGridParser& deck);
explicit GridManager(const Opm::EclipseGridParser& deck);
/// Construct a 3d corner-point grid or tensor grid from a deck.
explicit GridManager(Opm::DeckConstPtr newParserDeck);
/// Construct a 2d cartesian grid with cells of unit size.
GridManager(int nx, int ny);
@ -60,7 +65,7 @@ namespace Opm
/// Construct a grid from an input file.
/// The file format used is currently undocumented,
/// and is therefore only suited for internal use.
GridManager(const std::string& input_filename);
explicit GridManager(const std::string& input_filename);
/// Destructor.
~GridManager();
@ -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_;

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

@ -48,7 +48,7 @@ private:
/// Psuedo-constructor, can appear in template
template <typename Format> unique_ptr <OutputWriter>
create (const ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid) {
return unique_ptr <OutputWriter> (new Format (params, parser, grid));
}
@ -61,7 +61,7 @@ create (const ParameterGroup& params,
/// to the list below!
typedef map <const char*, unique_ptr <OutputWriter> (*)(
const ParameterGroup&,
std::shared_ptr <const EclipseGridParser>,
std::shared_ptr <EclipseGridParser>,
std::shared_ptr <const UnstructuredGrid>)> map_t;
map_t FORMATS = {
{ "output_ecl", &create <EclipseWriter> },
@ -71,7 +71,7 @@ map_t FORMATS = {
unique_ptr <OutputWriter>
OutputWriter::create (const ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid) {
// allocate a list which will be filled with writers. this list
// is initially empty (no output).

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

@ -108,7 +108,7 @@ public:
*/
static std::unique_ptr <OutputWriter>
create (const parameter::ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid);
};

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

@ -236,6 +236,8 @@ namespace Opm
filterIntegerField("ACTNUM", new_ACTNUM_);
filterDoubleField("PORO", new_PORO_);
filterDoubleField("NTG", new_NTG_);
filterDoubleField("SWCR", new_SWCR_);
filterDoubleField("SOWCR", new_SOWCR_);
filterDoubleField("PERMX", new_PERMX_);
filterDoubleField("PERMY", new_PERMY_);
filterDoubleField("PERMZ", new_PERMZ_);
@ -264,6 +266,8 @@ namespace Opm
if (!new_ACTNUM_.empty()) sp.setIntegerField("ACTNUM", new_ACTNUM_);
if (!new_PORO_.empty()) sp.setFloatingPointField("PORO", new_PORO_);
if (!new_NTG_.empty()) sp.setFloatingPointField("NTG", new_NTG_);
if (!new_SWCR_.empty()) sp.setFloatingPointField("SWCR", new_SWCR_);
if (!new_SOWCR_.empty()) sp.setFloatingPointField("SOWCR", new_SOWCR_);
if (!new_PERMX_.empty()) sp.setFloatingPointField("PERMX", new_PERMX_);
if (!new_PERMY_.empty()) sp.setFloatingPointField("PERMY", new_PERMY_);
if (!new_PERMZ_.empty()) sp.setFloatingPointField("PERMZ", new_PERMZ_);
@ -291,12 +295,16 @@ namespace Opm
outputField(out, new_ACTNUM_, "ACTNUM");
outputField(out, new_PORO_, "PORO");
if (hasNTG()) {outputField(out, new_NTG_, "NTG");}
if (hasSWCR()) {outputField(out, new_SWCR_, "SWCR");}
if (hasSOWCR()) {outputField(out, new_SOWCR_, "SOWCR");}
outputField(out, new_PERMX_, "PERMX");
outputField(out, new_PERMY_, "PERMY");
outputField(out, new_PERMZ_, "PERMZ");
outputField(out, new_SATNUM_, "SATNUM");
}
bool hasNTG() const {return !new_NTG_.empty(); }
bool hasSWCR() const {return !new_SWCR_.empty(); }
bool hasSOWCR() const {return !new_SOWCR_.empty(); }
private:
EclipseGridParser parser_;
@ -309,6 +317,8 @@ namespace Opm
std::vector<int> new_ACTNUM_;
std::vector<double> new_PORO_;
std::vector<double> new_NTG_;
std::vector<double> new_SWCR_;
std::vector<double> new_SOWCR_;
std::vector<double> new_PERMX_;
std::vector<double> new_PERMY_;
std::vector<double> new_PERMZ_;

1
opm/core/io/eclipse/EclipseGridParser.cpp
View File

@ -126,6 +126,7 @@ namespace EclipseKeywords
string("EQUIL"), string("PVCDO"), string("TSTEP"),
string("PLYVISC"), string("PLYROCK"), string("PLYADS"),
string("PLYMAX"), string("TLMIXPAR"), string("WPOLYMER"),
string("PLYSHEAR"),
string("GRUPTREE"), string("GCONINJE"), string("GCONPROD"),
string("WGRUPCON"), string("ENDSCALE"), string("SCALECRS"),
string("ENPTVD"), string("ENKRVD"),

1
opm/core/io/eclipse/EclipseGridParser.hpp
View File

@ -195,6 +195,7 @@ public:
SPECIAL_FIELD(PLYROCK)
SPECIAL_FIELD(PLYADS)
SPECIAL_FIELD(PLYMAX)
SPECIAL_FIELD(PLYSHEAR)
SPECIAL_FIELD(TLMIXPAR)
SPECIAL_FIELD(WPOLYMER)
SPECIAL_FIELD(GRUPTREE)

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

File diff suppressed because it is too large Load Diff

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

@ -24,7 +24,10 @@
#include <opm/core/io/OutputWriter.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <string>
#include <vector>
#include <memory> // std::unique_ptr
struct UnstructuredGrid;
@ -57,7 +60,15 @@ public:
* binary files using ERT.
*/
EclipseWriter(const parameter::ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid);
/*!
* \brief Sets the common attributes required to write eclipse
* binary files using ERT.
*/
EclipseWriter(const parameter::ParameterGroup& params,
Opm::DeckConstPtr newParserDeck,
std::shared_ptr <const UnstructuredGrid> grid);
/**
@ -87,16 +98,23 @@ public:
private:
std::shared_ptr <const EclipseGridParser> parser_;
Opm::DeckConstPtr newParserDeck_;
std::shared_ptr <const UnstructuredGrid> grid_;
bool enableOutput_;
int outputInterval_;
int outputTimeStepIdx_;
std::string outputDir_;
std::string baseName_;
PhaseUsage uses_; // active phases in the input deck
std::shared_ptr <EclipseSummary> summary_;
void activeToGlobalCellData_(std::vector<double> &globalCellsBuf,
const std::vector<double> &activeCellsBuf,
const std::vector<double> &inactiveCellsBuf) const;
/// Write solution field variables (pressure and saturation)
void writeSolution (const SimulatorTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState);
void writeSolution_(const SimulatorTimer& timer,
const SimulatorState& reservoirState);
};
} // namespace Opm

36
opm/core/io/eclipse/SpecialEclipseFields.hpp
View File

@ -2094,6 +2094,42 @@ struct PLYMAX : public SpecialBase
};
struct PLYSHEAR : public SpecialBase
{
std::vector<double> water_velocity_;
std::vector<double> vrf_;
virtual std::string name() const {return std::string("PLYSHEAR");}
virtual void read(std::istream& is)
{
// Note. This function assumes that NTSFUN = 1, and reads only one table.
std::vector<double> plyshear;
readVectorData(is, plyshear);
for (int i=0; i<(int)plyshear.size(); i+=2) {
water_velocity_.push_back(plyshear[i]);
vrf_.push_back(plyshear[i+1]);
}
}
virtual void write(std::ostream& os) const
{
os << name() << '\n';
for (int i=0; i<(int)water_velocity_.size(); ++i) {
os << water_velocity_[i] << " " << vrf_[i] << '\n';
}
os << '\n';
}
virtual void convertToSI(const EclipseUnits& units)
{
for (int i=0; i<(int)water_velocity_.size(); ++i) {
water_velocity_[i] *= units.length / units.time;
}
}
};
struct TLMIXPAR : public SpecialBase
{
std::vector<double> tlmixpar_;

13
opm/core/linalg/LinearSolverIstl.cpp
View File

@ -30,6 +30,7 @@
// TODO: clean up includes.
#include <dune/common/deprecated.hh>
#include <dune/common/version.hh>
#include <dune/istl/bvector.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
@ -39,6 +40,10 @@
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
#if DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
#include <dune/istl/paamg/fastamg.hh>
#endif
#include <stdexcept>
#include <iostream>
@ -61,7 +66,7 @@ namespace Opm
solveCG_AMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps);
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double prolongateFactor, int smoothsteps);
@ -173,7 +178,7 @@ namespace Opm
linsolver_prolongate_factor_, linsolver_smooth_steps_);
break;
case KAMG:
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
res = solveKAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_, linsolver_smooth_steps_);
#else
@ -181,7 +186,7 @@ namespace Opm
#endif
break;
case FastAMG:
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
res = solveFastAMG(A, x, b, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_);
#else
@ -311,7 +316,7 @@ namespace Opm
}
#ifdef HAS_DUNE_FAST_AMG
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
LinearSolverInterface::LinearSolverReport
solveKAMG(const Mat& A, Vector& x, Vector& b, double tolerance, int maxit, int verbosity,
double linsolver_prolongate_factor, int linsolver_smooth_steps)

4
opm/core/pressure/tpfa/cfs_tpfa_residual.c
View File

@ -869,7 +869,7 @@ assemble_completion_to_well(int i, int w, int c, int nc, int np,
W = wells->W;
ctrl = W->ctrls[ w ];
if (well_controls_get_current(ctrl) < 0) {
if (well_controls_well_is_shut( ctrl )) {
/* Interpreting a negative current control index to mean a shut well */
welleq_coeff_shut(np, h, &res, &w2c, &w2w);
}
@ -933,7 +933,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells *wells ,
for (w = i = 0; w < W->number_of_wells; w++) {
pw = wpress[ w ];
is_open = (well_controls_get_current(W->ctrls[w]) >= 0);
is_open = well_controls_well_is_open(W->ctrls[w]);
for (; i < W->well_connpos[w + 1]; i++, pmobp += np) {

2
opm/core/pressure/tpfa/ifs_tpfa.c
View File

@ -363,7 +363,7 @@ assemble_well_contrib(int nc ,
for (w = 0; w < W->number_of_wells; w++) {
ctrls = W->ctrls[ w ];
if (well_controls_get_current(ctrls) < 0) {
if (well_controls_well_is_shut(ctrls) ) {
/* Treat this well as a shut well, isolated from the domain. */

88
opm/core/props/BlackoilPropertiesFromDeck.cpp
View File

@ -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=*/0);
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, grid, /*numSamples=*/0);
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,
@ -58,8 +76,8 @@ namespace Opm
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 0);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (deck.hasField("ENDSCALE") && threephase_model != "simple") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
if (deck.hasField("ENDSCALE") && threephase_model != "gwseg") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
@ -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", 0);
pvt_.init(newParserDeck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 0);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "gwseg");
if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "gwseg") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'gwseg' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, grid, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
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<SatFuncStone2Nonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, grid, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, grid, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
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()
{
}

29
opm/core/props/BlackoilPropertiesFromDeck.hpp
View File

@ -27,6 +27,9 @@
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <memory>
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();

76
opm/core/props/phaseUsageFromDeck.hpp
View File

@ -24,10 +24,49 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
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

88
opm/core/props/pvt/BlackoilPvtProperties.cpp
View File

@ -32,6 +32,9 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvdcoTable.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
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,90 @@ 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");
if (phase_usage_.phase_used[Liquid]) {
densities_[phase_usage_.phase_pos[Liquid]]
= densityKeyword->getRecord(region_number_)->getItem("OIL")->getSIDouble(0);
}
if (phase_usage_.phase_used[Aqua]) {
densities_[phase_usage_.phase_pos[Aqua]]
= densityKeyword->getRecord(region_number_)->getItem("WATER")->getSIDouble(0);
}
if (phase_usage_.phase_used[Vapour]) {
densities_[phase_usage_.phase_pos[Vapour]]
= densityKeyword->getRecord(region_number_)->getItem("GAS")->getSIDouble(0);
}
} 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_);
if (samples > 0) {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
} else {
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDead(pvdoTable));
}
} 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_);
if (samples > 0) {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
} else {
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDead(pvdgTable));
}
} 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_;

10
opm/core/props/pvt/BlackoilPvtProperties.hpp
View File

@ -20,11 +20,12 @@
#ifndef OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
#define OPM_BLACKOILPVTPROPERTIES_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <string>
#include <memory>
@ -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;

31
opm/core/props/pvt/SinglePvtConstCompr.hpp
View File

@ -20,9 +20,12 @@
#ifndef OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
#define OPM_SINGLEPVTCONSTCOMPR_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/parser/eclipse/Utility/PvtwTable.hpp>
#include <opm/parser/eclipse/Utility/PvdcoTable.hpp>
#include <vector>
#include <algorithm>
@ -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),

48
opm/core/props/pvt/SinglePvtDead.cpp
View File

@ -62,6 +62,54 @@ namespace Opm
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
/// Constructor
SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
{
// Copy data
const std::vector<double>& press = pvdoTable.getPressureColumn();
const std::vector<double>& b = pvdoTable.getFormationFactorColumn();
const std::vector<double>& visc = pvdoTable.getViscosityColumn();
const int sz = b.size();
std::vector<double> bInv(sz);
for (int i = 0; i < sz; ++i) {
bInv[i] = 1.0 / b[i];
}
b_ = NonuniformTableLinear<double>(press, bInv);
viscosity_ = NonuniformTableLinear<double>(press, visc);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
/// Constructor
SinglePvtDead::SinglePvtDead(const Opm::PvdgTable& pvdgTable)
{
// Copy data
const std::vector<double>& press = pvdgTable.getPressureColumn();
const std::vector<double>& b = pvdgTable.getFormationFactorColumn();
const std::vector<double>& visc = pvdgTable.getViscosityColumn();
const int sz = b.size();
std::vector<double> bInv(sz);
for (int i = 0; i < sz; ++i) {
bInv[i] = 1.0 / b[i];
}
b_ = NonuniformTableLinear<double>(press, bInv);
viscosity_ = NonuniformTableLinear<double>(press, visc);
// Dumping the created tables.
// static int count = 0;
// std::ofstream os((std::string("dump-") + boost::lexical_cast<std::string>(count++)).c_str());
// os.precision(15);
// os << "1/B\n\n" << one_over_B_
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
// Destructor
SinglePvtDead::~SinglePvtDead()
{

6
opm/core/props/pvt/SinglePvtDead.hpp
View File

@ -23,6 +23,10 @@
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
#include <opm/parser/eclipse/Utility/PvdgTable.hpp>
#include <vector>
namespace Opm
@ -40,6 +44,8 @@ namespace Opm
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDead(const table_t& pvd_table);
SinglePvtDead(const Opm::PvdoTable& pvdoTable);
SinglePvtDead(const Opm::PvdgTable& pvdgTable);
virtual ~SinglePvtDead();
/// Viscosity as a function of p and z.

45
opm/core/props/pvt/SinglePvtDeadSpline.cpp
View File

@ -20,6 +20,9 @@
#include "config.h"
#include <opm/core/props/pvt/SinglePvtDeadSpline.hpp>
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/parser/eclipse/Utility/SimpleTable.hpp>
#include <algorithm>
// 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<std::string>(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<double> &press = pvdoTable.getPressureColumn();
const std::vector<double> &B = pvdoTable.getFormationFactorColumn();
const std::vector<double> &visc = pvdoTable.getViscosityColumn();
std::vector<double> 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<double> &press = pvdgTable.getPressureColumn();
const std::vector<double> &B = pvdgTable.getFormationFactorColumn();
const std::vector<double> &visc = pvdgTable.getViscosityColumn();
std::vector<double> 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

7
opm/core/props/pvt/SinglePvtDeadSpline.hpp
View File

@ -20,9 +20,12 @@
#ifndef OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#define OPM_SINGLEPVTDEADSPLINE_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
#include <opm/parser/eclipse/Utility/PvdgTable.hpp>
#include <vector>
namespace Opm
@ -41,6 +44,8 @@ namespace Opm
typedef std::vector<std::vector<std::vector<double> > > 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.

24
opm/core/props/pvt/SinglePvtLiveGas.cpp
View File

@ -32,8 +32,8 @@
#include <opm/core/props/pvt/SinglePvtLiveGas.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
#include <algorithm>
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; i<sz; ++i) {
const auto &undersatTable = *pvtgTable.getInnerTable(i);
undersat_gas_tables_[i].resize(3);
undersat_gas_tables_[i][0] = undersatTable.getOilSolubilityColumn();
undersat_gas_tables_[i][1] = undersatTable.getGasFormationFactorColumn();
undersat_gas_tables_[i][2] = pvtgTable.getOuterTable()->getGasViscosityColumn();
}
}
// Destructor
SinglePvtLiveGas::~SinglePvtLiveGas()
{

4
opm/core/props/pvt/SinglePvtLiveGas.hpp
View File

@ -21,6 +21,9 @@
#define OPM_SINGLEPVTLIVEGAS_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/parser/eclipse/Utility/PvtgTable.hpp>
#include <vector>
namespace Opm
@ -38,6 +41,7 @@ namespace Opm
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtLiveGas(const table_t& pvto);
SinglePvtLiveGas(const Opm::PvtgTable& pvtgTable);
virtual ~SinglePvtLiveGas();
/// Viscosity as a function of p and z.

69
opm/core/props/pvt/SinglePvtLiveOil.cpp
View File

@ -21,8 +21,8 @@
#include <opm/core/props/pvt/SinglePvtLiveOil.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
#include <algorithm>
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; i<sz; ++i) {
saturated_oil_table_[0][i] = saturatedPvto->getPressureColumn()[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; i<sz; ++i) {
const auto undersaturatedPvto = pvtoTable.getInnerTable(i);
undersat_oil_tables_[i].resize(3);
int tsize = undersaturatedPvto->numRows();
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; j<tsize; ++j) {
undersat_oil_tables_[i][0][j] = undersaturatedPvto->getPressureColumn()[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<sz; ++i) {
// Skip records already containing undersaturated data
if (undersat_oil_tables_[i][0].size() > 1) {
continue;
}
// Look ahead for next record containing undersaturated data
if (iNext < i) {
iNext = i+1;
while (iNext<sz && undersat_oil_tables_[iNext][0].size() < 2) {
++iNext;
}
if (iNext == sz) OPM_THROW(std::runtime_error,"Unable to complete undersaturated table.");
}
// Add undersaturated data to current record while maintaining compressibility and viscosibility
typedef std::vector<std::vector<std::vector<double> > >::size_type sz_t;
for (sz_t j=1; j<undersat_oil_tables_[iNext][0].size(); ++j) {
double diffPressure = undersat_oil_tables_[iNext][0][j]-undersat_oil_tables_[iNext][0][j-1];
double pressure = undersat_oil_tables_[i][0].back()+diffPressure;
undersat_oil_tables_[i][0].push_back(pressure);
double compr = (1.0/undersat_oil_tables_[iNext][1][j]-1.0/undersat_oil_tables_[iNext][1][j-1])
/ (0.5*(1.0/undersat_oil_tables_[iNext][1][j]+1.0/undersat_oil_tables_[iNext][1][j-1]));
double B = (1.0/undersat_oil_tables_[i][1].back())*(1.0+0.5*compr)/(1.0-0.5*compr);
undersat_oil_tables_[i][1].push_back(1.0/B);
double visc = (undersat_oil_tables_[iNext][2][j]-undersat_oil_tables_[iNext][2][j-1])
/ (0.5*(undersat_oil_tables_[iNext][2][j]+undersat_oil_tables_[iNext][2][j-1]));
double mu = (undersat_oil_tables_[i][2].back())*(1.0+0.5*visc)/(1.0-0.5*visc);
undersat_oil_tables_[i][2].push_back(mu);
}
}
}
/// Destructor.

5
opm/core/props/pvt/SinglePvtLiveOil.hpp
View File

@ -20,8 +20,10 @@
#ifndef OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
#define OPM_SINGLEPVTLIVEOIL_HEADER_INCLUDED
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/parser/eclipse/Utility/PvtoTable.hpp>
#include <vector>
namespace Opm
@ -39,6 +41,7 @@ namespace Opm
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtLiveOil(const table_t& pvto);
SinglePvtLiveOil(const Opm::PvtoTable &pvtoTable);
virtual ~SinglePvtLiveOil();
/// Viscosity as a function of p and z.

41
opm/core/props/rock/RockCompressibility.cpp
View File

@ -25,6 +25,9 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Utility/RocktabTable.hpp>
#include <opm/parser/eclipse/Utility/RockTable.hpp>
#include <iostream>
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);

6
opm/core/props/rock/RockCompressibility.hpp
View File

@ -20,6 +20,8 @@
#ifndef OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
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]

250
opm/core/props/rock/RockFromDeck.cpp
View File

@ -21,6 +21,9 @@
#include "config.h"
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <array>
namespace Opm
@ -37,6 +40,10 @@ namespace Opm
PermeabilityKind fillTensor(const EclipseGridParser& parser,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap);
PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& 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<double>& 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<const std::vector<double>*> tensor;
tensor.reserve(10);
const std::vector<double> zero(num_global_cells, 0.0);
tensor.push_back(&zero);
std::array<int,9> 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<unsigned char>::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<const std::vector<double>*>& tensor,
std::array<int,9>& 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

16
opm/core/props/rock/RockFromDeck.hpp
View File

@ -22,6 +22,9 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
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<double> porosity_;
std::vector<double> permeability_;

141
opm/core/props/satfunc/SatFuncBase.hpp
View File

@ -22,6 +22,11 @@
#include <opm/core/utility/buildUniformMonotoneTable.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Utility/SwofTable.hpp>
#include <opm/parser/eclipse/Utility/SgofTable.hpp>
#include <vector>
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 <class TableType>
void SatFuncBase<TableType>::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<double>& sw = swof.getSwColumn();
const std::vector<double>& krw = swof.getKrwColumn();
const std::vector<double>& krow = swof.getKrowColumn();
const std::vector<double>& 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<double> sw_ex(n+2);
std::vector<double> krw_ex(n+2);
std::vector<double> krow_ex(n+2);
std::vector<double> 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<double>::size_type i=1; i<sw.size(); ++i) {
if (krw[i]> 0.0) {
swcr_ = sw[i-1];
krorw_ = krow[i-1];
break;
}
}
for (std::vector<double>::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<double>& sg = sgof.getSgColumn();
const std::vector<double>& krg = sgof.getKrgColumn();
const std::vector<double>& krog = sgof.getKrogColumn();
const std::vector<double>& pcog = sgof.getPcogColumn();
// Extend the tables with constant values such that the
// derivatives at the endpoints are zero
int n = sg.size();
std::vector<double> sg_ex(n+2);
std::vector<double> krg_ex(n+2);
std::vector<double> krog_ex(n+2);
std::vector<double> 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<double>::size_type i=1; i<sg.size(); ++i) {
if (krg[i]> 0.0) {
sgcr_ = sg[i-1];
krorg_ = krog[i-1];
break;
}
}
for (std::vector<double>::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 <class TableType>
void SatFuncBase<TableType>::updateSatHyst(const double* s,
const EPSTransforms* epst,

10
opm/core/props/satfunc/SatFuncSimple.cpp
View File

@ -34,9 +34,9 @@ namespace Opm
template<>
void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int /*samples*/)
const std::vector<double>& arg,
const std::vector<double>& value,
const int /* samples */)
{
table = NonuniformTableLinear<double>(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;
@ -217,7 +217,7 @@ namespace Opm
out << "sg_shift: " << sg_shift << std::endl;
out << "sow_hyst: " << sow_hyst << std::endl;
out << "sow_shift: " << sow_shift << std::endl;
};
}
} // namespace Opm

16
opm/core/props/satfunc/SatFuncSimple.hpp
View File

@ -32,22 +32,22 @@ namespace Opm
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;
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<UniformTableLinear<double> > SatFuncSimpleUniform;
typedef SatFuncSimple<NonuniformTableLinear<double> > 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_); // ????

15
opm/core/props/satfunc/SatFuncStone2.hpp
View File

@ -32,24 +32,23 @@ namespace Opm
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<UniformTableLinear<double> > SatFuncStone2Uniform;
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;

27
opm/core/props/satfunc/SaturationPropsFromDeck.hpp
View File

@ -27,6 +27,9 @@
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
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<double>& 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<double>& scaleparam);
void initEPSParam(const int cell,
EPSTransforms::Transform& data,
const bool oil,
@ -149,6 +171,11 @@ namespace Opm
const std::vector<double>& s0,
const std::vector<double>& krsr,
const std::vector<double>& krmax);
bool columnIsMasked_(Opm::DeckConstPtr newParserDeck,
const std::string& keywordName,
int /* columnIdx */)
{ return newParserDeck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
};

464
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
View File

@ -26,6 +26,11 @@
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/parser/eclipse/Utility/EndscaleWrapper.hpp>
#include <opm/parser/eclipse/Utility/ScalecrsWrapper.hpp>
#include <opm/parser/eclipse/Utility/EnptvdTable.hpp>
#include <opm/parser/eclipse/Utility/EnkrvdTable.hpp>
#include <iostream>
namespace Opm
@ -140,6 +145,123 @@ namespace Opm
}
/// Initialize from deck.
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::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<int>& 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 <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid)
{
std::vector<double> swl, swcr, swu, sgl, sgcr, sgu, sowcr, sogcr;
std::vector<double> 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 <class SatFuncSet>
@ -450,6 +635,71 @@ namespace Opm
}
}
// Initialize hysteresis saturation scaling parameters
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSHyst(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid)
{
std::vector<double> iswl, iswcr, iswu, isgl, isgcr, isgu, isowcr, isogcr;
std::vector<double> 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 <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(const EclipseGridParser& deck,
@ -624,6 +874,213 @@ namespace Opm
}
}
// Initialize saturation scaling parameter
template <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSKey(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const std::string& keyword,
std::vector<double>& 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<std::vector<std::vector<double> > > table_dummy;
std::vector<std::vector<std::vector<double> > >& 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<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_aqua];
}
} else if (keyword == std::string("SWCR") || keyword == std::string("ISWCR") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 1))) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).swcr_;
}
} else if (keyword == std::string("SWU") || keyword == std::string("ISWU") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 2))) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_aqua];
}
} else if (keyword == std::string("SGL") || keyword == std::string("ISGL") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 3))) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smin_[phase_pos_vapour];
}
} else if (keyword == std::string("SGCR") || keyword == std::string("ISGCR") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 4))) {
itab = 5;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sgcr_;
}
} else if (keyword == std::string("SGU") || keyword == std::string("ISGU") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 5))) {
itab = 6;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).smax_[phase_pos_vapour];
}
} else if (keyword == std::string("SOWCR") || keyword == std::string("ISOWCR") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 6))) {
itab = 7;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sowcr_;
}
} else if (keyword == std::string("SOGCR") || keyword == std::string("ISOGCR") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENPTVD", 7))) {
itab = 8;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).sogcr_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 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<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwmax_;
}
} else if (keyword == std::string("KRG") || keyword == std::string("IKRG") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 1))) {
itab = 2;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgmax_;
}
if (useLiquid && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 2))) {
itab = 3;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).kromax_;
}
} else if (keyword == std::string("KRWR") || keyword == std::string("IKRWR") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 3))) {
itab = 4;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krwr_;
}
} else if (keyword == std::string("KRGR") || keyword == std::string("IKRGR") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 4))) {
itab = 5;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krgr_;
}
} else if (keyword == std::string("KRORW") || keyword == std::string("IKRORW") ) {
if (useAqua && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 5))) {
itab = 6;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorw_;
}
} else if (keyword == std::string("KRORG") || keyword == std::string("IKRORG") ) {
if (useVapour && (useKeyword || columnIsMasked_(newParserDeck, "ENKRVD", 6))) {
itab = 7;
scaleparam.resize(grid.number_of_cells);
for (int i=0; i<grid.number_of_cells; ++i)
scaleparam[i] = funcForCell(i).krorg_;
}
} else {
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 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<double>& 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<double>& depth = table[0][jtab];
std::vector<double>& 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 <class SatFuncSet>
void SaturationPropsFromDeck<SatFuncSet>::initEPSParam(const int cell,
@ -646,6 +1103,13 @@ namespace Opm
{
if (scr.empty() && su.empty() && (sxcr.empty() || !do_3pt_) && s0.empty()) {
data.doNotScale = true;
data.smin = sl_tab;
if (oil) {
data.smax = (s0_tab < 0.0) ? 1.0 - su_tab : 1.0 - su_tab - s0_tab;
} else {
data.smax = su_tab;
}
data.scr = scr_tab;
} else {
data.doNotScale = false;
data.do_3pt = do_3pt_;

8
opm/core/props/satfunc/SaturationPropsInterface.hpp
View File

@ -73,6 +73,14 @@ namespace Opm
const int* cells,
double* smin,
double* smax) const = 0;
/// Update saturation state for the hysteresis tracking
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
virtual void updateSatHyst(const int n,
const int* cells,
const double* s) = 0;
};

4
opm/core/simulator/SimulatorCompressibleTwophase.cpp
View File

@ -498,13 +498,13 @@ namespace Opm
<< std::endl;
std::cout.precision(8);
watercut.push(timer.currentTime() + timer.currentStepLength(),
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
if (wells_) {
wellreport.push(props_, *wells_,
state.pressure(), state.surfacevol(), state.saturation(),
timer.currentTime() + timer.currentStepLength(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
sreport.total_time = step_timer.secsSinceStart();

6
opm/core/simulator/SimulatorIncompTwophase.cpp
View File

@ -310,7 +310,7 @@ namespace Opm
const int nw = wells->number_of_wells;
for (int w = 0; w < nw; ++w) {
const WellControls* wc = wells->ctrls[w];
if (well_controls_get_current(wc) >= 0) {
if (well_controls_well_is_open( wc )) {
if (well_controls_get_current_type(wc) == BHP) {
return false;
}
@ -578,12 +578,12 @@ namespace Opm
dynamic_cast<TransportSolverTwophaseReorder&>(*tsolver_)
.solveGravity(&initial_porevol[0], stepsize, state);
}
watercut.push(timer.currentTime() + timer.currentStepLength(),
watercut.push(timer.simulationTimeElapsed() + timer.currentStepLength(),
produced[0]/(produced[0] + produced[1]),
tot_produced[0]/tot_porevol_init);
if (wells_) {
wellreport.push(props_, *wells_, state.saturation(),
timer.currentTime() + timer.currentStepLength(),
timer.simulationTimeElapsed() + timer.currentStepLength(),
well_state.bhp(), well_state.perfRates());
}
}

4
opm/core/simulator/SimulatorOutput.cpp
View File

@ -30,7 +30,7 @@ using namespace Opm;
SimulatorOutputBase::SimulatorOutputBase (
const parameter::ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid,
std::shared_ptr <const SimulatorTimer> timer,
std::shared_ptr <const SimulatorState> state,
@ -72,7 +72,7 @@ SimulatorOutputBase::operator std::function <void ()> () {
void
SimulatorOutputBase::writeOutput () {
const int this_time = timer_->currentTime ();
const int this_time = timer_->simulationTimeElapsed ();
// if the simulator signals for timesteps that aren't reporting
// times, then ignore them

6
opm/core/simulator/SimulatorOutput.hpp
View File

@ -53,7 +53,7 @@ protected:
* need to pick them up from the object members.
*/
SimulatorOutputBase (const parameter::ParameterGroup& p,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid,
std::shared_ptr <const SimulatorTimer> timer,
std::shared_ptr <const SimulatorState> state,
@ -139,7 +139,7 @@ private:
template <typename Simulator>
struct SimulatorOutput : public SimulatorOutputBase {
SimulatorOutput (const parameter::ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid,
std::shared_ptr <const SimulatorTimer> timer,
std::shared_ptr <const SimulatorState> state,
@ -161,7 +161,7 @@ struct SimulatorOutput : public SimulatorOutputBase {
* the arguments passed exceeds the lifetime of this object.
*/
SimulatorOutput (const parameter::ParameterGroup& params,
const EclipseGridParser& parser,
EclipseGridParser& parser,
const UnstructuredGrid& grid,
const SimulatorTimer& timer,
const SimulatorState& state,

100
opm/core/simulator/SimulatorTimer.cpp
View File

@ -57,28 +57,65 @@ namespace Opm
start_date_ = deck.getStartDate();
}
/// Use the SimulatorTimer as a shim around opm-parser's Opm::TimeMap
void SimulatorTimer::init(Opm::TimeMapConstPtr timeMap,
size_t beginReportStepIdx,
size_t endReportStepIdx)
{
timeMap_ = timeMap;
current_step_ = 0;
beginReportStepIdx_ = beginReportStepIdx;
endReportStepIdx_ = std::min(timeMap_->numTimesteps() + 1, endReportStepIdx);
}
/// Total number of steps.
int SimulatorTimer::numSteps() const
{
return timesteps_.size();
if (timeMap_)
return endReportStepIdx_ - beginReportStepIdx_;
else
return timesteps_.size();
}
/// Index of the first considered simulation episode
size_t SimulatorTimer::beginReportStepIndex() const
{
if (!timeMap_) {
OPM_THROW(std::runtime_error, "indexFirstEpisode() is only implemented "
"for simulation timers which are based on Opm::TimeMap");
}
return beginReportStepIdx_;
}
/// Index of the last considered simulation episode
size_t SimulatorTimer::endReportStepIndex() const
{
if (!timeMap_) {
OPM_THROW(std::runtime_error, "indexLastEpisode() is only implemented "
"for simulation timers which are based on Opm::TimeMap");
}
return endReportStepIdx_;
}
/// Current step number.
int SimulatorTimer::currentStepNum() const
{
return current_step_;
return current_step_ + beginReportStepIdx_;
}
/// 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 +123,44 @@ namespace Opm
double SimulatorTimer::currentStepLength() const
{
assert(!done());
return timesteps_[current_step_];
if (timeMap_)
return timeMap_->getTimeStepLength(beginReportStepIdx_ + 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(beginReportStepIdx_ + current_step_ - 1);
else
return timesteps_[current_step_ - 1];
}
/// Current time.
double SimulatorTimer::currentTime() const
/// time elapsed since the start of the simulation [s].
double SimulatorTimer::simulationTimeElapsed() const
{
return current_time_;
if (timeMap_)
return
timeMap_->getTimePassedUntil(beginReportStepIdx_ + current_step_);
else
return current_time_;
}
/// time elapsed since the start of the POSIX epoch (Jan 1st, 1970) [s].
time_t SimulatorTimer::currentPosixTime() const
{
tm t = boost::posix_time::to_tm(currentDateTime());
return std::mktime(&t);
}
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(beginReportStepIdx_ + current_step_);
else
return boost::posix_time::ptime(start_date_) + boost::posix_time::seconds( (int) current_time_ );
}
@ -112,7 +168,11 @@ namespace Opm
/// Total time.
double SimulatorTimer::totalTime() const
{
return total_time_;
if (timeMap_)
return
timeMap_->getTotalTime();
else
return total_time_;
}
/// Set total time.
@ -121,14 +181,20 @@ 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.
void SimulatorTimer::report(std::ostream& os) const
{
os << "\n\n--------------- Simulation step number " << currentStepNum() << " ---------------"
<< "\n Current time (days) " << Opm::unit::convert::to(currentTime(), Opm::unit::day)
<< "\n Current time (days) " << Opm::unit::convert::to(simulationTimeElapsed(), Opm::unit::day)
<< "\n Current stepsize (days) " << Opm::unit::convert::to(currentStepLength(), Opm::unit::day)
<< "\n Total time (days) " << Opm::unit::convert::to(totalTime(), Opm::unit::day)
<< "\n" << std::endl;
@ -138,7 +204,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 +213,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 current_step_ > int(endReportStepIdx_ - beginReportStepIdx_ - 1);
else
return int(timesteps_.size()) == current_step_;
}

25
opm/core/simulator/SimulatorTimer.hpp
View File

@ -20,6 +20,8 @@
#ifndef OPM_SIMULATORTIMER_HEADER_INCLUDED
#define OPM_SIMULATORTIMER_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp>
#include <iosfwd>
#include <vector>
#include <boost/date_time/gregorian/gregorian.hpp>
@ -47,9 +49,20 @@ 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,
size_t beginReportStepIdx = 0,
size_t endReportStepIdx = std::numeric_limits<size_t>::max());
/// Total number of steps.
int numSteps() const;
/// Index of the first report step considered
size_t beginReportStepIndex() const;
/// Index of the next-after-last report step to be considered
size_t endReportStepIndex() const;
/// Current step number. This is the number of timesteps that
/// has been completed from the start of the run. The time
/// after initialization but before the simulation has started
@ -73,8 +86,13 @@ namespace Opm
/// it is an error to call stepLengthTaken().
double stepLengthTaken () const;
/// Current time.
double currentTime() const;
/// Time elapsed since the start of the POSIX epoch (Jan 1st,
/// 1970) until the current time step begins [s].
time_t currentPosixTime() const;
/// Time elapsed since the start of the simulation until the
/// beginning of the current time step [s].
double simulationTimeElapsed() const;
/// Return the current time as a posix time object.
boost::posix_time::ptime currentDateTime() const;
@ -99,7 +117,10 @@ namespace Opm
bool done() const;
private:
Opm::TimeMapConstPtr timeMap_;
std::vector<double> timesteps_;
size_t beginReportStepIdx_;
size_t endReportStepIdx_;
int current_step_;
double current_time_;
double total_time_;

68
opm/core/simulator/WellState.hpp
View File

@ -23,6 +23,7 @@
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
#include <vector>
#include <cassert>
namespace Opm
{
@ -45,28 +46,55 @@ namespace Opm
bhp_.resize(nw);
wellrates_.resize(nw * np, 0.0);
for (int w = 0; w < nw; ++w) {
assert((wells->type[w] == INJECTOR) || (wells->type[w] == PRODUCER));
const WellControls* ctrl = wells->ctrls[w];
// Initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a little
// above or below (depending on if the well is an
// injector or producer) pressure in first perforation
// cell.
if ((well_controls_get_current(ctrl) < 0) || // SHUT
(well_controls_get_current_type(ctrl) != BHP)) {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*state.pressure()[first_cell];
} else {
bhp_[w] = well_controls_get_current_target( ctrl );
}
// Initialize well rates to match controls if type is SURFACE_RATE
if ((well_controls_get_current(ctrl) >= 0) && // open well
(well_controls_get_current_type(ctrl) == SURFACE_RATE)) {
const double rate_target = well_controls_get_current_target(ctrl);
const double * distr = well_controls_get_current_distr( ctrl );
if (well_controls_well_is_shut(ctrl)) {
// Shut well:
// 1. Assign zero well rates.
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = rate_target * distr[p];
wellrates_[np*w + p] = 0.0;
}
// 2. Assign bhp equal to bhp control, if
// applicable, otherwise assign equal to
// first perforation cell pressure.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp_[w] = well_controls_get_current_target( ctrl );
} else {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
bhp_[w] = state.pressure()[first_cell];
}
} else {
// Open well:
// 1. Initialize well rates to match controls
// if type is SURFACE_RATE. Otherwise, we
// cannot set the correct value here, so we
// assign a small rate with the correct
// sign so that any logic depending on that
// sign will work as expected.
if (well_controls_get_current_type(ctrl) == SURFACE_RATE) {
const double rate_target = well_controls_get_current_target(ctrl);
const double * distr = well_controls_get_current_distr( ctrl );
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = rate_target * distr[p];
}
} else {
const double small_rate = 1e-14;
const double sign = (wells->type[w] == INJECTOR) ? 1.0 : -1.0;
for (int p = 0; p < np; ++p) {
wellrates_[np*w + p] = small_rate * sign;
}
}
// 2. Initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
if (well_controls_get_current_type(ctrl) == BHP) {
bhp_[w] = well_controls_get_current_target( ctrl );
} else {
const int first_cell = wells->well_cells[wells->well_connpos[w]];
const double safety_factor = (wells->type[w] == INJECTOR) ? 1.01 : 0.99;
bhp_[w] = safety_factor*state.pressure()[first_cell];
}
}
}

145
opm/core/simulator/initState_impl.hpp
View File

@ -32,6 +32,8 @@
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
#include <iostream>
#include <cmath>
@ -471,6 +473,108 @@ namespace Opm
}
/// Initialize a state from input deck.
template <class Props, class State>
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") && newParserDeck->hasKeyword("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
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<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& 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<double>& 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<double>& 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<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
const std::vector<double>& 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 <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
@ -485,6 +589,10 @@ namespace Opm
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
if (deck.hasField("EQUIL") && deck.hasField("PRESSURE")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): The deck must either specify the initial "
"condition using the PRESSURE _or_ the EQUIL keyword (currently it has both)");
}
state.init(grid, num_phases);
if (deck.hasField("EQUIL")) {
if (num_phases != 2) {
@ -568,10 +676,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 +874,39 @@ namespace Opm
}
}
/// Initialize a blackoil state from input deck.
template <class Props, class State>
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<double>& 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<double>& 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

18
opm/core/well_controls.h
View File

@ -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);
@ -55,8 +55,18 @@ well_controls_get_current( const struct WellControls * ctrl);
void
well_controls_set_current( struct WellControls * ctrl, int current);
void
well_controls_invert_current( struct WellControls * ctrl );
bool
well_controls_well_is_shut(const struct WellControls * ctrl);
bool
well_controls_well_is_open(const struct WellControls * ctrl);
void
well_controls_open_well( struct WellControls * ctrl);
void
well_controls_shut_well( struct WellControls * ctrl);
int
well_controls_add_new(enum WellControlType type , double target , const double * distr , struct WellControls * ctrl);

2
opm/core/wells.h
View File

@ -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);

59
opm/core/wells/WellCollection.cpp
View File

@ -19,10 +19,66 @@
#include "config.h"
#include <opm/core/wells/WellCollection.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <boost/lexical_cast.hpp>
#include <memory>
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<WellsGroupInterface> child = createGroupWellsGroup(groupChild, timeStep, phaseUsage);
WellsGroup* parent_as_group = static_cast<WellsGroup*> (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<std::string>(timeStep) << " to group named " << wellChild->getGroupName(timeStep) << ", but this group does not exist in the WellCollection.");
}
std::shared_ptr<WellsGroupInterface> child = createWellWellsGroup(wellChild, timeStep, phaseUsage);
WellsGroup* parent_as_group = static_cast<WellsGroup*> (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<WellNode*>(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<WellsGroupInterface> 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<WellNode*>& WellCollection::getLeafNodes() const {
return leaf_nodes_;
}

15
opm/core/wells/WellCollection.hpp
View File

@ -23,10 +23,15 @@
#define OPM_WELLCOLLECTION_HPP
#include <vector>
#include <memory>
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/core/grid.h>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <memory>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
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.

60
opm/core/wells/WellsGroup.cpp
View File

@ -74,7 +74,7 @@ namespace Opm
{
return parent_;
}
const std::string& WellsGroupInterface::name()
const std::string& WellsGroupInterface::name() const
{
return name_;
}
@ -747,9 +747,7 @@ namespace Opm
void WellNode::shutWell()
{
if (shut_well_) {
// We set the tilde of the current control
// set_current_control(self_index_, -1, wells_);
well_controls_invert_current(wells_->ctrls[self_index_]);
well_controls_shut_well( wells_->ctrls[self_index_]);
}
else {
const double target = 0.0;
@ -767,7 +765,7 @@ namespace Opm
well_controls_iset_target( wells_->ctrls[self_index_] , group_control_index_ , target);
well_controls_iset_distr(wells_->ctrls[self_index_] , group_control_index_ , distr);
}
well_controls_invert_current(wells_->ctrls[self_index_]);
well_controls_open_well( wells_->ctrls[self_index_]);
}
}
@ -1141,4 +1139,56 @@ namespace Opm
return return_value;
}
std::shared_ptr<WellsGroupInterface> 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<WellsGroupInterface> wells_group(new WellsGroup(group->name(), production_specification, injection_specification, phase_usage));
return wells_group;
}
std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep, const PhaseUsage& phase_usage )
{
InjectionSpecification injection_specification;
ProductionSpecification production_specification;
if (well->isInjector(timeStep)) {
const WellInjectionProperties& properties = well->getInjectionProperties(timeStep);
injection_specification.BHP_limit_ = properties.BHPLimit;
injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(properties.injectorType));
injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(properties.controlMode));
injection_specification.surface_flow_max_rate_ = properties.surfaceInjectionRate;
injection_specification.reservoir_flow_max_rate_ = properties.reservoirInjectionRate;
production_specification.guide_rate_ = 0.0; // We know we're not a producer
}
else if (well->isProducer(timeStep)) {
const WellProductionProperties& properties = well->getProductionProperties(timeStep);
production_specification.BHP_limit_ = properties.BHPLimit;
production_specification.reservoir_flow_max_rate_ = properties.ResVRate;
production_specification.oil_max_rate_ = properties.OilRate;
production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(properties.controlMode));
production_specification.water_max_rate_ = properties.WaterRate;
injection_specification.guide_rate_ = 0.0; // we know we're not an injector
}
std::shared_ptr<WellsGroupInterface> wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage));
return wells_group;
}
}

20
opm/core/wells/WellsGroup.hpp
View File

@ -25,6 +25,10 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <string>
#include <memory>
@ -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<WellsGroupInterface> 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<WellsGroupInterface> 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<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep,
const PhaseUsage& phase_usage );
}
#endif /* OPM_WELLSGROUP_HPP */

1036
opm/core/wells/WellsManager.cpp
View File

File diff suppressed because it is too large Load Diff

92
opm/core/wells/WellsManager.hpp
View File

@ -21,8 +21,11 @@
#define OPM_WELLSMANAGER_HEADER_INCLUDED
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/core/wells/WellCollection.hpp>
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
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,35 @@ 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.
explicit 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::EclipseStateConstPtr eclipseState,
const size_t timeStep,
const UnstructuredGrid& grid,
const double* permeability);
/// Destructor.
~WellsManager();
/// 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 +127,31 @@ namespace Opm
void applyExplicitReinjectionControls(const std::vector<double>& well_reservoirrates_phase,
const std::vector<double>& 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<int,int>& cartesian_to_compressed );
void setupWellControls(std::vector<WellConstPtr>& wells, size_t timeStep,
std::vector<std::string>& well_names, const PhaseUsage& phaseUsage);
// Data
Wells* w_;
void createWellsFromSpecs( std::vector<WellConstPtr>& wells, size_t timeStep,
const UnstructuredGrid& grid,
std::vector<std::string>& well_names,
std::vector<WellData>& well_data,
std::map<std::string, int> & well_names_to_index,
const PhaseUsage& phaseUsage,
const std::map<int,int> cartesian_to_compressed,
const double* permeability);
void addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage);
void setupGuideRates(std::vector<WellConstPtr>& wells, const size_t timeStep, std::vector<WellData>& well_data, std::map<std::string, int>& well_names_to_index);
// Data
Wells* w_;
WellCollection well_collection_;
};
} // namespace Opm

54
opm/core/wells/well_controls.c
View File

@ -98,6 +98,8 @@ struct WellControls
*/
int current;
bool well_is_open;
/*
The capacity allocated.
*/
@ -130,7 +132,7 @@ well_controls_create(void)
ctrl = malloc(1 * sizeof *ctrl);
if (ctrl != NULL) {
/* Initialise empty control set */
/* Initialise empty control set; the well is created open. */
ctrl->num = 0;
ctrl->number_of_phases = 0;
ctrl->type = NULL;
@ -138,6 +140,7 @@ well_controls_create(void)
ctrl->distr = NULL;
ctrl->current = -1;
ctrl->cpty = 0;
ctrl->well_is_open = true;
}
return ctrl;
@ -192,11 +195,23 @@ well_controls_set_current( struct WellControls * ctrl, int current) {
ctrl->current = current;
}
void
well_controls_invert_current( struct WellControls * ctrl ) {
ctrl->current = ~ctrl->current;
bool well_controls_well_is_shut(const struct WellControls * ctrl) {
return !ctrl->well_is_open;
}
bool well_controls_well_is_open(const struct WellControls * ctrl) {
return ctrl->well_is_open;
}
void well_controls_open_well( struct WellControls * ctrl) {
ctrl->well_is_open = true;
}
void well_controls_shut_well( struct WellControls * ctrl) {
ctrl->well_is_open = false;
}
enum WellControlType
well_controls_iget_type(const struct WellControls * ctrl, int control_index) {
@ -293,19 +308,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;

37
opm/core/wells/wells.c
View File

@ -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");
}
{

4
tests/testBlackoilState1.DATA
View File

@ -10,5 +10,9 @@ DYV
DZV
10*10 /
ACTNUM
1 998*2 3 /
DEPTHZ
121*2000 /

186
tests/test_linearsolver.cpp Normal file
View File

@ -0,0 +1,186 @@
/*
Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE OPM-IterativeSolverTest
#include <boost/test/unit_test.hpp>
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <dune/common/version.hh>
#include <memory>
#include <cstdlib>
#include <string>
struct MyMatrix
{
MyMatrix(int rows, int nnz)
: data(nnz, 0.0), rowStart(rows+1, -1),
colIndex(nnz, -1)
{}
MyMatrix()
: data(), rowStart(), colIndex()
{}
std::vector<double> data;
std::vector<int> rowStart;
std::vector<int> colIndex;
};
std::shared_ptr<MyMatrix> createLaplacian(int N)
{
MyMatrix* mm=new MyMatrix(N*N, N*N*5);
int nnz=0;
mm->rowStart[0]=0;
for(int row=0; row<N*N; row++)
{
int x=row%N;
int y=row/N;
if(y>0)
{
mm->colIndex[nnz]=row-N;
mm->data[nnz++]=-1;
}
if(x>0)
{
mm->colIndex[nnz]=row-1;
mm->data[nnz++]=-1;
}
mm->colIndex[nnz]=row;
mm->data[nnz++]=4;
if(x<N-1)
{
mm->colIndex[nnz]=row+1;
mm->data[nnz++]=-1;
}
if(y<N-1)
{
mm->colIndex[nnz]=row+N;
mm->data[nnz++]=-1;
}
mm->rowStart[row+1]=nnz;
}
mm->data.resize(nnz);
mm->colIndex.resize(nnz);
return std::shared_ptr<MyMatrix>(mm);
}
void createRandomVectors(int NN, std::vector<double>& x, std::vector<double>& b,
const MyMatrix& mat)
{
x.resize(NN);
for(auto entry=x.begin(), end =x.end(); entry!=end; ++entry)
*entry=((double) (rand()%100))/10.0;
b.resize(NN);
std::fill(b.begin(), b.end(), 0.0);
// Construct the right hand side as b=A*x
for(std::size_t row=0; row<mat.rowStart.size()-1; ++row)
{
for(int i=mat.rowStart[row], end=mat.rowStart[row+1]; i!=end; ++i)
{
b[row]+= mat.data[i]*x[mat.colIndex[i]];
}
}
}
void run_test(const Opm::parameter::ParameterGroup& param)
{
int N=4;
auto mat = createLaplacian(N);
std::vector<double> x(N*N), b(N*N);
createRandomVectors(100*100, x, b, *mat);
std::vector<double> exact(x);
std::fill(x.begin(), x.end(), 0.0);
Opm::LinearSolverFactory ls(param);
ls.solve(N*N, mat->data.size(), &(mat->rowStart[0]),
&(mat->colIndex[0]), &(mat->data[0]), &(b[0]),
&(x[0]));
}
BOOST_AUTO_TEST_CASE(DefaultTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
#ifdef HAVE_DUNE_ISTL
BOOST_AUTO_TEST_CASE(CGAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("1"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(CGILUTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("0"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(BiCGILUTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("2"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
BOOST_AUTO_TEST_CASE(FastAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("3"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
param.insertParameter(std::string("linsolver_verbosity"), std::string("2"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(KAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("4"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
run_test(param);
}
#endif
#endif

54
tests/test_parser.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <boost/test/unit_test.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/Deck/DeckDoubleItem.hpp>
#include <string>
#include <iostream>
#include <vector>
#include <memory>
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<int>& 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] );
}

91
tests/test_wellcollection.cpp Normal file
View File

@ -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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <boost/test/unit_test.hpp>
#include <opm/core/wells/WellCollection.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTreeNode.hpp>
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<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
for (size_t i=0; i<wells.size(); i++) {
collection.addWell(wells[i], 2, pu);
}
BOOST_CHECK_EQUAL("G1", collection.findNode("INJ1")->getParent()->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());
}

23
tests/test_wellcontrols.cpp
View File

@ -45,11 +45,6 @@ BOOST_AUTO_TEST_CASE(Construction)
well_controls_set_current( ctrls , 2 );
BOOST_CHECK_EQUAL( 2 , well_controls_get_current( ctrls ));
well_controls_invert_current( ctrls );
BOOST_CHECK( well_controls_get_current( ctrls ) < 0 );
well_controls_invert_current( ctrls );
BOOST_CHECK_EQUAL( 2 , well_controls_get_current( ctrls ));
{
enum WellControlType type1 = BHP;
enum WellControlType type2 = SURFACE_RATE;
@ -103,3 +98,21 @@ BOOST_AUTO_TEST_CASE(Construction)
}
BOOST_AUTO_TEST_CASE(OpenClose)
{
struct WellControls * ctrls = well_controls_create();
BOOST_CHECK_EQUAL( true , well_controls_well_is_open(ctrls) );
BOOST_CHECK_EQUAL( false , well_controls_well_is_shut(ctrls) );
well_controls_open_well( ctrls );
BOOST_CHECK_EQUAL( true , well_controls_well_is_open(ctrls) );
BOOST_CHECK_EQUAL( false , well_controls_well_is_shut(ctrls) );
well_controls_shut_well( ctrls );
BOOST_CHECK_EQUAL( false , well_controls_well_is_open(ctrls) );
BOOST_CHECK_EQUAL( true , well_controls_well_is_shut(ctrls) );
well_controls_destroy( ctrls );
}

6
tests/test_wells.cpp
View File

@ -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<Wells> 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<Wells> 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 ));
}

110
tests/test_wellsgroup.cpp Normal file
View File

@ -0,0 +1,110 @@
/*
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 <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#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 <memory>
#include <vector>
#include <boost/test/unit_test.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTreeNode.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
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<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
for (size_t i=0; i<wells.size(); i++) {
WellConstPtr well = wells[i];
std::shared_ptr<WellsGroupInterface> wellsGroup = createWellWellsGroup(well, 2, pu);
BOOST_CHECK_EQUAL(well->name(), wellsGroup->name());
if (well->isInjector(2)) {
const WellInjectionProperties& properties = well->getInjectionProperties(2);
BOOST_CHECK_EQUAL(properties.surfaceInjectionRate, wellsGroup->injSpec().surface_flow_max_rate_);
BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->injSpec().BHP_limit_);
BOOST_CHECK_EQUAL(properties.reservoirInjectionRate, wellsGroup->injSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_);
}
if (well->isProducer(2)) {
const WellProductionProperties& properties = well->getProductionProperties(2);
BOOST_CHECK_EQUAL(properties.ResVRate, wellsGroup->prodSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(properties.BHPLimit, wellsGroup->prodSpec().BHP_limit_);
BOOST_CHECK_EQUAL(properties.OilRate, wellsGroup->prodSpec().oil_max_rate_);
BOOST_CHECK_EQUAL(properties.WaterRate, 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<GroupTreeNodeConstPtr> nodes = eclipseState->getSchedule()->getGroupTree(2)->getNodes();
for (size_t i=0; i<nodes.size(); i++) {
GroupConstPtr group = eclipseState->getSchedule()->getGroup(nodes[i]->name());
std::shared_ptr<WellsGroupInterface> 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_);
}
}
}

93
tests/test_wellsmanager.cpp
View File

@ -27,6 +27,10 @@
#define BOOST_TEST_MODULE WellsManagerTests
#include <boost/test/unit_test.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Schedule.hpp>
#include <opm/core/wells/WellsManager.hpp>
#include <opm/core/wells.h>
#include <opm/core/well_controls.h>
@ -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
@ -169,29 +173,24 @@ void check_controls_epoch1( struct WellControls ** 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")));
BOOST_AUTO_TEST_CASE(Constructor_Works) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
const Wells* wells = wellsManager.c_wells();
wells_static_check( wells );
check_controls_epoch0( wells->ctrls );
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
wells_static_check( wellsManager.c_wells() );
check_controls_epoch0( wellsManager.c_wells()->ctrls );
}
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(Deck, *gridManager.c_grid(), NULL);
const Wells* wells = wellsManager.c_wells();
wells_static_check( wells );
check_controls_epoch1( wells->ctrls );
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
wells_static_check( wellsManager.c_wells() );
check_controls_epoch1( wellsManager.c_wells()->ctrls );
}
}
@ -200,15 +199,14 @@ BOOST_AUTO_TEST_CASE(Constructor_Works) {
BOOST_AUTO_TEST_CASE(WellsEqual) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
Deck.setCurrentEpoch(0);
Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager1(eclipseState , 1 , *gridManager.c_grid(), NULL);
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) );
}
@ -216,21 +214,36 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
Opm::EclipseGridParser Deck("wells_manager_data.data");
Opm::GridManager gridManager(Deck);
Deck.setCurrentEpoch(0);
Opm::WellsManager wellsManager0(Deck, *gridManager.c_grid(), NULL);
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data.data")));
Deck.setCurrentEpoch(1);
Opm::WellsManager wellsManager1(Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager0(eclipseState , 0 , *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager1(eclipseState , 1 , *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 );
}

68
tests/wells_group.data Executable file
View File

@ -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

30
tests/wells_manager_data.data
View File

@ -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

86
tests/wells_manager_data_expanded.data Executable file
View File

@ -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

42
tests/wells_manager_data_wellSTOP.data Executable file
View File

@ -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