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

Merge remote-tracking branch 'origin/opm-parser-integrate' into refactor-for-cpgrid-support

Browse Source 
Resolved Conflicts:
	opm/core/props/BlackoilPropertiesFromDeck.cpp
	opm/core/props/rock/RockFromDeck.hpp
	opm/core/props/satfunc/SaturationPropsFromDeck.hpp
	opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
This commit is contained in:
Markus Blatt
2014-02-19 15:22:18 +01:00
parent 8da99e11c5 f3cbaafb65
commit 3aad1a6cc5
61 changed files with 4489 additions and 2009 deletions
Download Patch File Download Diff File Expand all files Collapse all files

25
CMakeLists_files.cmake
View File

@@ -166,11 +166,13 @@ list (APPEND TEST_SOURCE_FILES
tests/test_param.cpp
tests/test_blackoilfluid.cpp
tests/test_shadow.cpp
tests/test_units.cpp
tests/test_blackoilstate.cpp
tests/test_parser.cpp
tests/test_wellsmanager.cpp
tests/test_wellcontrols.cpp
tests/test_units.cpp
tests/test_blackoilstate.cpp
tests/test_parser.cpp
tests/test_wellsmanager.cpp
tests/test_wellcontrols.cpp
tests/test_wellsgroup.cpp
tests/test_wellcollection.cpp
)
# originally generated with the command:
@@ -179,11 +181,13 @@ list (APPEND TEST_DATA_FILES
tests/extratestdata.xml
tests/testdata.xml
tests/liveoil.DATA
tests/wetgas.DATA
tests/testBlackoilState1.DATA
tests/testBlackoilState2.DATA
tests/wells_manager_data.data
tests/wells_manager_data_expanded.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:
@@ -309,6 +313,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/rock/RockBasic.hpp
opm/core/props/rock/RockCompressibility.hpp
opm/core/props/rock/RockFromDeck.hpp
opm/core/props/satfunc/SatFuncBase.hpp
opm/core/props/satfunc/SatFuncGwseg.hpp
opm/core/props/satfunc/SatFuncSimple.hpp
opm/core/props/satfunc/SatFuncStone2.hpp

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"

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

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

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

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

@@ -19,5 +19,5 @@ set (opm-autodiff_DEPS
dune-istl REQUIRED;
opm-core REQUIRED"
# Eigen
"Eigen3 3.1 REQUIRED"
"Eigen3 3.1.2 "
)

3
examples/compute_tof.cpp
View File

@@ -68,7 +68,6 @@ namespace
}
void buildTracerheadsFromWells(const Wells* wells,
const Opm::WellState& well_state,
Opm::SparseTable<int>& tracerheads)
{
if (wells == 0) {
@@ -257,7 +256,7 @@ try
std::vector<double> tracer;
Opm::SparseTable<int> tracerheads;
if (compute_tracer) {
buildTracerheadsFromWells(wells->c_wells(), well_state, tracerheads);
buildTracerheadsFromWells(wells->c_wells(), tracerheads);
}
if (use_dg) {
if (compute_tracer) {

167
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
{
@@ -216,6 +302,76 @@ namespace Opm
}
top_depths_vec.resize(x.size()*y.size(), tops[0]);
top_depths = &top_depths_vec[0];
} else {
OPM_THROW(std::runtime_error,
"Could not find either TOPS or DEPTHZ keyword, "
"one of them is required for initialization with DXV/DYV/DZV.");
}
// 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.");
}
}
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.
@@ -227,5 +383,4 @@ namespace Opm
}
} // namespace Opm

11
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
{
@@ -44,6 +46,9 @@ namespace Opm
/// Construct a 3d corner-point grid or tensor grid from a deck.
GridManager(const Opm::EclipseGridParser& deck);
/// Construct a 3d corner-point grid or tensor grid from a deck.
GridManager(Opm::DeckConstPtr newParserDeck);
/// Construct a 2d cartesian grid with cells of unit size.
GridManager(int nx, int ny);
@@ -72,6 +77,8 @@ namespace Opm
/// to make it clear that we are returning a C-compatible struct.
const UnstructuredGrid* c_grid() const;
static void createGrdecl(Opm::DeckConstPtr newParserDeck, struct grdecl &grdecl);
private:
// Disable copying and assignment.
GridManager(const GridManager& other);
@@ -79,8 +86,10 @@ namespace Opm
// Construct corner-point grid from deck.
void initFromDeckCornerpoint(const Opm::EclipseGridParser& deck);
void initFromDeckCornerpoint(Opm::DeckConstPtr newParserDeck);
// Construct tensor grid from deck.
void initFromDeckTensorgrid(const Opm::EclipseGridParser& deck);
void initFromDeckTensorgrid(Opm::DeckConstPtr newParserDeck);
// The managed UnstructuredGrid.
UnstructuredGrid* ug_;

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

@@ -100,10 +100,17 @@ namespace EclipseKeywords
string("RV"),
string("DXV"), string("DYV"), string("DZV"),
string("DEPTHZ"), string("TOPS"), string("MAPAXES"),
string("SWCR"), string("SWL"), string("SWU"),
string("SOWCR"), string("KRW"), string("KRWR"),
string("KRO"), string("KRORW"), string("NTG"),
string("RHO")
string("SWL"), string("SWCR"), string("SWU"),
string("SGL"), string("SGCR"), string("SGU"),
string("SOWCR"), string("SOGCR"), string("KRW"),
string("KRWR"), string("KRG"), string("KRGR"),
string("KRO"), string("KRORW"), string("KRORG"),
string("ISWL"), string("ISWCR"), string("ISWU"),
string("ISGL"), string("ISGCR"), string("ISGU"),
string("ISOWCR"), string("ISOGCR"), string("IKRW"),
string("IKRWR"), string("IKRG"), string("IKRGR"),
string("IKRO"), string("IKRORW"), string("IKRORG"),
string("NTG"), string("RHO")
};
const int num_floating_fields = sizeof(floating_fields) / sizeof(floating_fields[0]);
@@ -561,9 +568,17 @@ void EclipseGridParser::convertToSI()
key == "LAMEMOD" || key == "SHEARMOD" || key == "POISSONMOD" ||
key == "PWAVEMOD" || key == "MULTPV" || key == "PWAVEMOD" ||
key == "SGAS" || key == "SWAT" || key == "SOIL" ||
key == "NTG" || key == "SWCR" || key == "SWL" ||
key == "SWU" || key == "SOWCR" || key == "KRW" ||
key == "KRWR" || key == "KRORW" || key == "KRO" ||
key == "SWL" || key == "SWCR" || key == "SWU" ||
key == "SGL" || key == "SGCR" || key == "SGU" ||
key == "SOWCR" || key == "SOGCR" || key == "KRW" ||
key == "KRWR" || key == "KRG" || key == "KRGR" ||
key == "KRO" || key == "KRORW" || key == "KRORG" ||
key == "ISWL" || key == "ISWCR" || key == "ISWU" ||
key == "ISGL" || key == "ISGCR" || key == "ISGU" ||
key == "ISOWCR" || key == "ISOGCR" || key == "IKRW" ||
key == "IKRWR" || key == "IKRG" || key == "IKRGR" ||
key == "IKRO" || key == "IKRORW" || key == "IKRORG" ||
key == "NTG" ||
key == "RHO") /* nonstandard field with no unit logic. use with caution */ {
unit = 1.0;
do_convert = false; // Dimensionless keywords...

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

@@ -575,6 +575,18 @@ private:
fortio_fclose) { }
};
// in order to get RTTI for this "class" (which is just a typedef), we must
// ask the compiler to explicitly instantiate it.
template struct EclipseHandle<ecl_sum_tstep_struct>;
} // anonymous namespace
// Note: the following parts were taken out of the anonymous
// namespace, since EclipseSummary is now used as a pointer member in
// EclipseWriter and forward declared in EclipseWriter.hpp.
// forward decl. of mutually dependent type
struct EclipseWellReport;
@@ -655,9 +667,6 @@ private:
}
};
// in order to get RTTI for this "class" (which is just a typedef), we must
// ask the compiler to explicitly instantiate it.
template struct EclipseHandle<ecl_sum_tstep_struct>;
/**
* Summary variable that reports a characteristics of a well.
@@ -670,7 +679,7 @@ protected:
PhaseUsage uses, /* phases present */
BlackoilPhases::PhaseIndex phase, /* oil, water or gas */
WellType type, /* prod. or inj. */
char aggregation, /* rate or total */
char aggregation, /* rate or total or BHP */
std::string unit)
: EclipseHandle <smspec_node_type> (
ecl_sum_add_var (summary,
@@ -718,22 +727,26 @@ private:
char aggregation) {
std::string name;
name += 'W'; // well
switch (phase) {
if (aggregation == 'B') {
name += "BHP";
} else {
switch (phase) {
case BlackoilPhases::Aqua: name += 'W'; break; /* water */
case BlackoilPhases::Vapour: name += 'G'; break; /* gas */
case BlackoilPhases::Liquid: name += 'O'; break; /* oil */
default:
OPM_THROW(std::runtime_error,
"Unknown phase used in blackoil reporting");
}
switch (type) {
}
switch (type) {
case WellType::INJECTOR: name += 'I'; break;
case WellType::PRODUCER: name += 'P'; break;
default:
OPM_THROW(std::runtime_error,
"Unknown well type used in blackoil reporting");
}
name += aggregation; /* rate ('R') or total ('T') */
}
name += aggregation; /* rate ('R') or total ('T') */
return name;
}
protected:
@@ -743,6 +756,13 @@ protected:
const double value = sign_ * wellState.wellRates () [index_] * convFactor;
return value;
}
double bhp (const WellState& wellstate) {
// Note that 'index_' is used here even though it is meant
// to give a (well,phase) pair.
const int num_phases = wellstate.wellRates().size() / wellstate.bhp().size();
return wellstate.bhp()[index_/num_phases];
}
};
/// Monitors the rate given by a well.
@@ -761,7 +781,7 @@ struct EclipseWellRate : public EclipseWellReport {
type,
'R',
"SM3/DAY" /* surf. cub. m. per day */ ) { }
virtual double update (const SimulatorTimer& timer,
virtual double update (const SimulatorTimer& /*timer*/,
const WellState& wellState) {
// TODO: Why only positive rates?
return std::max (0., rate (wellState));
@@ -790,6 +810,11 @@ struct EclipseWellTotal : public EclipseWellReport {
virtual double update (const SimulatorTimer& timer,
const WellState& wellState) {
if (timer.currentStepNum() == 0) {
// We are at the initial state.
// No step has been taken yet.
return 0.0;
}
// TODO: Is the rate average for the timestep, or is in
// instantaneous (in which case trapezoidal or Simpson integration
// would probably be better)
@@ -805,6 +830,31 @@ private:
double total_;
};
/// Monitors the bottom hole pressure in a well.
struct EclipseWellBhp : public EclipseWellReport {
EclipseWellBhp (const EclipseSummary& summary,
const EclipseGridParser& parser,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
parser,
whichWell,
uses,
phase,
type,
'B',
"Pascal")
{ }
virtual double update (const SimulatorTimer& /*timer*/,
const WellState& wellState)
{
return bhp(wellState);
}
};
inline void
EclipseSummary::writeTimeStep (const SimulatorTimer& timer,
const WellState& wellState) {
@@ -861,8 +911,30 @@ EclipseSummary::addWells (const EclipseGridParser& parser,
}
}
}
// Add BHP monitors
for (int whichWell = 0; whichWell != numWells; ++whichWell) {
// In the call below: uses, phase and the well type arguments
// are not used, except to set up an index that stores the
// well indirectly. For details see the implementation of the
// EclipseWellReport constructor, and the method
// EclipseWellReport::bhp().
BlackoilPhases::PhaseIndex phase = BlackoilPhases::Liquid;
if (!uses.phase_used[BlackoilPhases::Liquid]) {
phase = BlackoilPhases::Vapour;
}
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellBhp (*this,
parser,
whichWell,
uses,
phase,
WELL_TYPES[0])));
}
}
namespace {
/// Helper method that can be used in keyword transformation (must curry
/// the barsa argument)
static double toBar (const double& pressure) {
@@ -902,11 +974,16 @@ void EclipseWriter::writeInit(const SimulatorTimer &timer,
/* Initial solution (pressure and saturation) */
writeSolution (timer, reservoirState, wellState);
/* Create summary object (could not do it at construction time,
since it requires knowledge of the start time). */
summary_.reset(new EclipseSummary(outputDir_, baseName_, timer, *parser_));
summary_->addWells (*parser_, uses_);
}
void EclipseWriter::writeSolution (const SimulatorTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState) {
const WellState& /*wellState*/) {
// start writing to files
EclipseRestart rst (outputDir_,
baseName_,
@@ -952,9 +1029,15 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// (first timestep, in practice), and reused later. but how to do this
// without keeping the complete summary in memory (which will then
// accumulate all the timesteps)?
EclipseSummary sum (outputDir_, baseName_, timer, *parser_);
sum.addWells (*parser_, uses_);
sum.writeTimeStep (timer, wellState);
//
// Note: The answer to the question above is still not settled,
// but now we do keep the complete summary in memory, as a member
// variable in the EclipseWriter class, instead of creating a
// temporary EclipseSummary in this function every time it is
// called. This has been changed so that the final summary file
// will contain data from the whole simulation, instead of just
// the last step.
summary_->writeTimeStep (timer, wellState);
}
#else

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

@@ -28,6 +28,7 @@
#include <memory> // std::unique_ptr
struct UnstructuredGrid;
struct EclipseSummary;
namespace Opm {
@@ -90,6 +91,7 @@ private:
std::string outputDir_;
std::string baseName_;
PhaseUsage uses_; // active phases in the input deck
std::shared_ptr <EclipseSummary> summary_;
/// Write solution field variables (pressure and saturation)
void writeSolution (const SimulatorTimer& timer,

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

@@ -2313,8 +2313,8 @@ struct ENPTVD : public SpecialBase {
}
virtual void read(std::istream & is) {
table_.resize(5);
std::vector<std::vector<double> > sub_table(5);
table_.resize(9);
std::vector<std::vector<double> > sub_table(9);
while (!is.eof()) {
if(is.peek() == int('/')) {
if (sub_table[0].empty() && !(table_[0].empty())) {
@@ -2332,19 +2332,14 @@ struct ENPTVD : public SpecialBase {
if (data[0] == -1.0) {
OPM_THROW(std::runtime_error, "Error reading ENPTVD data - depth can not be defaulted.");
}
if ((data[4] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0) || (data[8] != -1.0)) {
OPM_THROW(std::runtime_error, "Error reading ENPTVD data - non-default values in column 5-7,9 not supported.");
for(std::vector<std::vector<double> >::size_type i=0; i<sub_table.size(); ++i) {
sub_table[i].push_back(data[i]); // [0-8]: depth swl swcr swu sgl sgcr sgu sowcr sogcr
}
sub_table[0].push_back(data[0]); //depth
sub_table[1].push_back(data[1]); //swl
sub_table[2].push_back(data[2]); //swcr
sub_table[3].push_back(data[3]); //swu
sub_table[4].push_back(data[7]); //sowcr
is >> ignoreWhitespace;
if(is.peek() == int('/')) {
is >> ignoreLine;
if (sub_table[0].size() >= 2) {
insertDefaultValues(sub_table, 5, -1.0, false);
insertDefaultValues(sub_table, 9, -1.0, false);
std::vector<std::vector<double> >::iterator it_sub = sub_table.begin();
for(std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
table_[i].push_back(*it_sub);
@@ -2369,9 +2364,9 @@ struct ENPTVD : public SpecialBase {
virtual void write(std::ostream & os) const {
os << name() << '\n';
std::cout << "-----depth-------swl------swcr-------swu-----sowcr" << std::endl;
std::cout << "-----depth-------swl------swcr-------swu-------sgl------sgcr-------sgu-----sowcr-----sogcr" << std::endl;
for (std::vector<std::vector<double> >::size_type j=0; j<table_[0].size(); ++j) {
std::cout << "--------------------------------------------------" << std::endl;
std::cout << "------------------------------------------------------------------------------------------" << std::endl;
for (std::vector<double>::size_type k=0; k<table_[0][j].size(); ++k) {
for (std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
std::cout << std::setw(10) << table_[i][j][k];
@@ -2400,8 +2395,8 @@ struct ENKRVD : public SpecialBase {
}
virtual void read(std::istream & is) {
table_.resize(5);
std::vector<std::vector<double> > sub_table(5);
table_.resize(8);
std::vector<std::vector<double> > sub_table(8);
while (!is.eof()) {
if(is.peek() == int('/')) {
if (sub_table[0].empty() && !(table_[0].empty())) {
@@ -2419,19 +2414,15 @@ struct ENKRVD : public SpecialBase {
if (data[0] == -1.0) {
OPM_THROW(std::runtime_error, "Error reading ENKRVD data - depth can not be defaulted.");
}
if ((data[2] != -1.0) || (data[5] != -1.0) || (data[6] != -1.0)) {
OPM_THROW(std::runtime_error, "Error reading ENKRVD data - non-default values in column 3,6-7 not supported.");
for(std::vector<std::vector<double> >::size_type i=0; i<sub_table.size(); ++i) {
sub_table[i].push_back(data[i]); // [0-7]: depth krw krg kro krwr krgr krorw krorg
}
sub_table[0].push_back(data[0]); //depth
sub_table[1].push_back(data[1]); //krw
sub_table[2].push_back(data[3]); //kro
sub_table[3].push_back(data[4]); //krw(sowcr)
sub_table[4].push_back(data[7]); //kro(swcr)
is >> ignoreWhitespace;
if(is.peek() == int('/')) {
is >> ignoreLine;
if (sub_table[0].size() >= 2) {
insertDefaultValues(sub_table, 5, -1.0, false);
insertDefaultValues(sub_table, 8, -1.0, false);
std::vector<std::vector<double> >::iterator it_sub = sub_table.begin();
for(std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
table_[i].push_back(*it_sub);
@@ -2457,9 +2448,9 @@ struct ENKRVD : public SpecialBase {
virtual void write(std::ostream & os) const {
os << name() << '\n';
std::cout << "-----depth-------krw------krow------krwr-----krorw" << std::endl;
std::cout << "-----depth-------krw-------krg-------kro------krwr------krgr-----krorw-----krorg" << std::endl;
for (std::vector<std::vector<double> >::size_type j=0; j<table_[0].size(); ++j) {
std::cout << "--------------------------------------------------" << std::endl;
std::cout << "--------------------------------------------------------------------------------" << std::endl;
for (std::vector<double>::size_type k=0; k<table_[0][j].size(); ++k) {
for (std::vector<std::vector<std::vector<double> > >::size_type i=0; i<table_.size(); ++i) {
std::cout << std::setw(10) << table_[i][j][k];

2
opm/core/linalg/LinearSolverUmfpack.cpp
View File

@@ -56,7 +56,7 @@ namespace Opm
const_cast<double*>(sa)
};
call_UMFPACK(&A, rhs, solution);
LinearSolverReport rep = {0};
LinearSolverReport rep = {};
rep.converged = true;
return rep;
}

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

@@ -23,7 +23,6 @@
namespace Opm
{
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
bool init_rock)
@@ -33,6 +32,14 @@ namespace Opm
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
bool init_rock)
{
init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims,
grid.cell_centroids, grid.dimensions, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
@@ -42,6 +49,15 @@ namespace Opm
grid.dimensions, param, init_rock);
}
BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock)
{
init(newParserDeck, grid.number_of_cells, grid.global_cell, grid.cartdims, grid.cell_centroids,
grid.dimensions, param, init_rock);
}
BlackoilPropertiesFromDeck::~BlackoilPropertiesFromDeck()
{
}

46
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();
@@ -202,6 +229,23 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock);
RockFromDeck rock_;
BlackoilPvtProperties pvt_;
std::unique_ptr<SaturationPropsInterface> satprops_;

102
opm/core/props/BlackoilPropertiesFromDeck_impl.hpp
View File

@@ -24,6 +24,31 @@ namespace Opm
}
}
template<class T>
void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock)
{
if (init_rock){
rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
}
pvt_.init(newParserDeck, /*numSamples=*/200);
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids, dimension,
/*numSamples=*/200);
if (pvt_.numPhases() != satprops_->numPhases()) {
OPM_THROW(std::runtime_error, "BlackoilPropertiesFromDeck::BlackoilPropertiesFromDeck() - Inconsistent number of phases in pvt data ("
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
template<class T>
void BlackoilPropertiesFromDeck::init(const EclipseGridParser& deck,
int number_of_cells,
@@ -98,4 +123,81 @@ namespace Opm
<< pvt_.numPhases() << ") and saturation-dependent function data (" << satprops_->numPhases() << ").");
}
}
template<class T>
void BlackoilPropertiesFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock)
{
if(init_rock){
rock_.init(newParserDeck, number_of_cells, global_cell, cart_dims);
}
const int pvt_samples = param.getDefault("pvt_tab_size", 200);
pvt_.init(newParserDeck, pvt_samples);
// Unfortunate lack of pointer smartness here...
const int sat_samples = param.getDefault("sat_tab_size", 200);
std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
if (newParserDeck->hasKeyword("ENDSCALE") && threephase_model != "simple") {
OPM_THROW(std::runtime_error, "Sorry, end point scaling currently available for the 'simple' model only.");
}
if (sat_samples > 1) {
if (threephase_model == "stone2") {
SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
= new SaturationPropsFromDeck<SatFuncStone2Uniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleUniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegUniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
dimension, 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, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "simple") {
SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
dimension, sat_samples);
} else if (threephase_model == "gwseg") {
SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
= new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
satprops_.reset(ptr);
ptr->init(newParserDeck, number_of_cells, global_cell, begin_cell_centroids,
dimension, 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() << ").");
}
}
}

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

@@ -24,13 +24,13 @@
#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)
@@ -105,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

81
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,83 @@ namespace Opm
}
}
void BlackoilPvtProperties::init(Opm::DeckConstPtr newParserDeck, int samples)
{
// If we need multiple regions, this class and the SinglePvt* classes must change.
region_number_ = 0;
phase_usage_ = phaseUsageFromDeck(newParserDeck);
// Surface densities. Accounting for different orders in eclipse and our code.
if (newParserDeck->hasKeyword("DENSITY")) {
Opm::DeckKeywordConstPtr densityKeyword = newParserDeck->getKeyword("DENSITY");
const std::vector<double>& d = densityKeyword->getRecord(region_number_)->getItem(0)->getSIDoubleData();
enum { ECL_oil = 0, ECL_water = 1, ECL_gas = 2 };
if (phase_usage_.phase_used[Aqua]) {
densities_[phase_usage_.phase_pos[Aqua]] = d[ECL_water];
}
if (phase_usage_.phase_used[Vapour]) {
densities_[phase_usage_.phase_pos[Vapour]] = d[ECL_gas];
}
if (phase_usage_.phase_used[Liquid]) {
densities_[phase_usage_.phase_pos[Liquid]] = d[ECL_oil];
}
} else {
OPM_THROW(std::runtime_error, "Input is missing DENSITY\n");
}
// Set the properties.
props_.resize(phase_usage_.num_phases);
// Water PVT
if (phase_usage_.phase_used[Aqua]) {
if (newParserDeck->hasKeyword("PVTW")) {
Opm::PvtwTable pvtwTable(newParserDeck->getKeyword("PVTW"));
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(pvtwTable));
} else {
// Eclipse 100 default.
props_[phase_usage_.phase_pos[Aqua]].reset(new SinglePvtConstCompr(0.5*Opm::prefix::centi*Opm::unit::Poise));
}
}
// Oil PVT
if (phase_usage_.phase_used[Liquid]) {
if (newParserDeck->hasKeyword("PVDO")) {
Opm::PvdoTable pvdoTable(newParserDeck->getKeyword("PVDO"), region_number_);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtDeadSpline(pvdoTable, samples));
} else if (newParserDeck->hasKeyword("PVTO")) {
Opm::PvtoTable pvtoTable(newParserDeck->getKeyword("PVTO"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtLiveOil(pvtoTable));
} else if (newParserDeck->hasKeyword("PVCDO")) {
Opm::PvdcoTable pvcdoTable(newParserDeck->getKeyword("PVCDO"));
props_[phase_usage_.phase_pos[Liquid]].reset(new SinglePvtConstCompr(pvcdoTable));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDO or PVTO\n");
}
}
// Gas PVT
if (phase_usage_.phase_used[Vapour]) {
if (newParserDeck->hasKeyword("PVDG")) {
Opm::PvdgTable pvdgTable(newParserDeck->getKeyword("PVDG"), region_number_);
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtDeadSpline(pvdgTable, samples));
} else if (newParserDeck->hasKeyword("PVTG")) {
Opm::PvtgTable pvtgTable(newParserDeck->getKeyword("PVTG"), /*tableIdx=*/0);
props_[phase_usage_.phase_pos[Vapour]].reset(new SinglePvtLiveGas(pvtgTable));
} else {
OPM_THROW(std::runtime_error, "Input is missing PVDG or PVTG\n");
}
}
// Must inform pvt property objects of phase structure.
for (int i = 0; i < phase_usage_.num_phases; ++i) {
props_[i]->setPhaseConfiguration(phase_usage_.num_phases, phase_usage_.phase_pos);
}
}
const double* BlackoilPvtProperties::surfaceDensities() const
{
return densities_;

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),

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

@@ -62,6 +62,30 @@ namespace Opm
// << "\n\nvisc\n\n" << viscosity_ << std::endl;
}
/// Constructor
SinglePvtDead::SinglePvtDead(const Opm::PvdoTable& pvdoTable)
{
// Copy data
const std::vector<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;
}
// Destructor
SinglePvtDead::~SinglePvtDead()
{

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

@@ -23,6 +23,9 @@
#include <opm/core/props/pvt/SinglePvtInterface.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/parser/eclipse/Utility/PvdoTable.hpp>
#include <vector>
namespace Opm
@@ -40,6 +43,7 @@ namespace Opm
public:
typedef std::vector<std::vector<std::vector<double> > > table_t;
SinglePvtDead(const table_t& pvd_table);
SinglePvtDead(const Opm::PvdoTable &pvdoTable);
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.

48
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()
{
@@ -99,12 +121,12 @@ namespace Opm
}
/// Viscosity and its derivatives as a function of p and r.
void SinglePvtLiveGas::mu(const int n,
const double* p,
const double* r,
double* output_mu,
double* output_dmudp,
double* output_dmudr) const
void SinglePvtLiveGas::mu(const int /*n*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_mu*/,
double* /*output_dmudp*/,
double* /*output_dmudr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");
}
@@ -156,12 +178,12 @@ namespace Opm
}
/// The inverse of the formation volume factor b = 1 / B, and its derivatives as a function of p and r.
void SinglePvtLiveGas::b(const int n,
const double* p,
const double* r,
double* output_b,
double* output_dbdp,
double* output_dbdr) const
void SinglePvtLiveGas::b(const int /*n*/,
const double* /*p*/,
const double* /*r*/,
double* /*output_b*/,
double* /*output_dbdp*/,
double* /*output_dbdr*/) const
{
OPM_THROW(std::runtime_error, "The new fluid interface not yet implemented");

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]

253
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
@@ -77,6 +84,17 @@ namespace Opm
assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold);
}
void RockFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell,
const int* cart_dims)
{
assignPorosity(newParserDeck, number_of_cells, global_cell);
permfield_valid_.assign(number_of_cells, false);
const double perm_threshold = 0.0; // Maybe turn into parameter?
assignPermeability(newParserDeck, number_of_cells, global_cell, cart_dims,
perm_threshold);
}
void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
int number_of_cells, const int* global_cell)
@@ -91,6 +109,20 @@ namespace Opm
}
}
void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell)
{
porosity_.assign(number_of_cells, 1.0);
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 = (global_cell == NULL) ? c : global_cell[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,
@@ -147,6 +179,61 @@ namespace Opm
}
}
void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cartdims,
double perm_threshold)
{
const int dim = 3;
const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
const int nc = 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 = 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 {
@@ -216,6 +303,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
@@ -345,6 +498,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

22
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;
@@ -52,7 +55,16 @@ namespace Opm
void init(const EclipseGridParser& deck,
int number_of_cells, const int* global_cell,
const int* cart_dims);
/// Initialize from deck and cell mapping.
/// \param newParserDeck Deck produced by the opm-parser code
/// \param number_of_cells The number of cells in the grid.
/// \param global_cell The mapping fom local to global cell indices.
/// global_cell[i] is the corresponding global index of i.
/// \param cart_dims The size of the underlying cartesian grid.
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell,
const int* cart_dims);
/// \return D, the number of spatial dimensions. Always 3 for deck input.
int numDimensions() const
{
@@ -83,11 +95,19 @@ namespace Opm
void assignPorosity(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell);
void assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell);
void assignPermeability(const EclipseGridParser& parser,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
const double perm_threshold);
void assignPermeability(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
double perm_threshold);
std::vector<double> porosity_;
std::vector<double> permeability_;

453
opm/core/props/satfunc/SatFuncBase.hpp Normal file
View File

@@ -0,0 +1,453 @@
/*
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef SATFUNCBASE_HPP
#define SATFUNCBASE_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#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
{
// Transforms for saturation table scaling
struct EPSTransforms {
struct Transform {
bool doNotScale;
bool do_3pt;
double smin;
double scr;
double sr;
double smax;
double slope1;
double slope2;
double scaleSat(double ss, double s_r, double s_cr, double s_max) const;
double scaleSatInv(double s, double s_r, double s_cr, double s_max) const;
double scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const; // Returns scaleSat'(s)
double scaleSatPc(double s, double s_min, double s_max) const;
double scaleSatDerivPc(double s, double s_min, double s_max) const; // Returns scaleSatPc'(s)
bool doKrMax;
bool doKrCrit;
bool doSatInterp;
double krsr;
double krmax;
double krSlopeMax;
double krSlopeCrit;
double scaleKr(double s, double kr, double krsr_tab) const;
double scaleKrDeriv(double s, double krDeriv) const; // Returns scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))
void printMe(std::ostream & out);
};
Transform wat;
Transform watoil;
Transform gas;
Transform gasoil;
};
// Hysteresis
struct SatHyst {
double sg_hyst;
double sg_shift;
double sow_hyst;
double sow_shift;
void printMe(std::ostream & out);
};
template <class TableType>
class SatFuncBase : public BlackoilPhases
{
public:
void init(const EclipseGridParser& deck,
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,
SatHyst* sat_hyst) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double krgmax_; // Max gas relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double sgcr_; // Critical gas saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double krgr_; // Gas relperm at critical oil-in-gas saturation.
double sowcr_; // Critical oil-in-water saturation.
double sogcr_; // Critical oil-in-gas-and-connate-water saturation.
double krorw_; // Oil relperm at critical water saturation.
double krorg_; // Oil relperm at critical gas saturation.
protected:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
TableType krw_;
TableType krow_;
TableType pcow_;
TableType krg_;
TableType krog_;
TableType pcog_;
double krocw_; // = krow_(s_wc)
private:
void extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
void initializeTableType(TableType& table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples);
};
template <class TableType>
void SatFuncBase<TableType>::init(const EclipseGridParser& deck,
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]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
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]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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>::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,
const EPSTransforms* epst_hyst,
SatHyst* sat_hyst) const
{
if (phase_usage.phase_used[Aqua] && phase_usage.phase_used[Vapour]) { //Water/Oil/Gas
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_-smin_[gpos], sowcr_, 1.0-smin_[wpos]-smin_[gpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_-smin_[wpos], sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
} else if (phase_usage.phase_used[Aqua]) { //Water/oil
int opos = phase_usage.phase_pos[BlackoilPhases::Liquid];
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
if (s[opos] > sat_hyst->sow_hyst)
{
sat_hyst->sow_hyst = s[opos];
double _sow_hyst = epst->watoil.scaleSat(sat_hyst->sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
double sow_hyst_shifted = epst_hyst->watoil.scaleSatInv(_sow_hyst, 1.0-swcr_, sowcr_, 1.0-smin_[wpos]);
sat_hyst->sow_shift = sow_hyst_shifted - sat_hyst->sow_hyst;
}
} else if (phase_usage.phase_used[Vapour]) {//Gas/Oil
int gpos = phase_usage.phase_pos[BlackoilPhases::Vapour];
if (s[gpos] > sat_hyst->sg_hyst)
{
sat_hyst->sg_hyst = s[gpos];
double _sg_hyst = epst->gas.scaleSat(sat_hyst->sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
double sg_hyst_shifted = epst_hyst->gas.scaleSatInv(_sg_hyst, 1.0-sogcr_, sgcr_, smax_[gpos]);
sat_hyst->sg_shift = sg_hyst_shifted - sat_hyst->sg_hyst;
}
}
}
template <class TableType>
void SatFuncBase<TableType>::extendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
} // namespace Opm
#endif // SATFUNCBASE_HPP

503
opm/core/props/satfunc/SatFuncGwseg.cpp
View File

@@ -29,504 +29,5 @@
namespace Opm
{
void SatFuncGwsegUniform::init(const EclipseGridParser& deck,
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]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
// 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);
SatFuncGwsegUniform::ExtendTable(sw,sw_ex,1);
SatFuncGwsegUniform::ExtendTable(krw,krw_ex,0);
SatFuncGwsegUniform::ExtendTable(krow,krow_ex,0);
SatFuncGwsegUniform::ExtendTable(pcow,pcow_ex,0);
buildUniformMonotoneTable(sw_ex, krw_ex, samples, krw_);
buildUniformMonotoneTable(sw_ex, krow_ex, samples, krow_);
buildUniformMonotoneTable(sw_ex, pcow_ex, samples, pcow_);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncGwsegUniform::ExtendTable(sg,sg_ex,1);
SatFuncGwsegUniform::ExtendTable(krg,krg_ex,0);
SatFuncGwsegUniform::ExtendTable(krog,krog_ex,0);
SatFuncGwsegUniform::ExtendTable(pcog,pcog_ex,0);
buildUniformMonotoneTable(sg_ex, krg_ex, samples, krg_);
buildUniformMonotoneTable(sg_ex, krog_ex, samples, krog_);
buildUniformMonotoneTable(sg_ex, pcog_ex, samples, pcog_);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncGwsegUniform::ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
void SatFuncGwsegUniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double xw = (sw - swco) / std::max(sg + sw - swco, eps);
const double xg = 1 - xw;
const double ssw = sg + sw;
const double ssg = sw - swco + sg;
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncGwsegUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
const double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
// xw and xg are the fractions occupied by water and gas zones.
const double eps = 1e-6;
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double krw = krw_(ssw);
const double krg = krg_(ssg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
const double xw = (sw - swco) / ssg;
const double xg = 1 - xw;
kr[Aqua] = xw*krw;
kr[Vapour] = xg*krg;
kr[Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = krw_.derivative(ssw);
const double dkrgg = krg_.derivative(ssg);
const double dkrow = krow_.derivative(ssw);
const double dkrog = krog_.derivative(ssg);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
dkrds[Aqua + Aqua*np] = (xg/d)*krw + xw*dkrww;
dkrds[Aqua + Vapour*np] = -(xw/d)*krw + xw*dkrww;
dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[Vapour + Aqua*np] = -(xg/d)*krg + xg*dkrgg;
dkrds[Vapour + Vapour*np] = (xw/d)*krg + xg*dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncGwsegUniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncGwsegUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncGwsegNonuniform ======
void SatFuncGwsegNonuniform::init(const EclipseGridParser& deck,
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]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
// 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);
SatFuncGwsegNonuniform::ExtendTable(sw,sw_ex,1);
SatFuncGwsegNonuniform::ExtendTable(krw,krw_ex,0);
SatFuncGwsegNonuniform::ExtendTable(krow,krow_ex,0);
SatFuncGwsegNonuniform::ExtendTable(pcow,pcow_ex,0);
krw_ = NonuniformTableLinear<double>(sw_ex, krw_ex);
krow_ = NonuniformTableLinear<double>(sw_ex, krow_ex);
pcow_ = NonuniformTableLinear<double>(sw_ex, pcow_ex);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncGwsegNonuniform::ExtendTable(sg,sg_ex,1);
SatFuncGwsegNonuniform::ExtendTable(krg,krg_ex,0);
SatFuncGwsegNonuniform::ExtendTable(krog,krog_ex,0);
SatFuncGwsegNonuniform::ExtendTable(pcog,pcog_ex,0);
krg_ = NonuniformTableLinear<double>(sg_ex, krg_ex);
krog_ = NonuniformTableLinear<double>(sg_ex, krog_ex);
pcog_ = NonuniformTableLinear<double>(sg_ex, pcog_ex);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncGwsegNonuniform::ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
void SatFuncGwsegNonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
// xw and xg are the fractions occupied by water and gas zones.
const double ssg = sw - swco + sg;
const double xw = (sw - swco) / ssg;
const double xg = 1 - xw;
const double ssw = sg + sw;
const double krw = krw_(sw);
const double krg = krg_(sg);
const double krow = krow_(ssw);
const double krog = krog_(ssg);
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncGwsegNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
double swco = smin_[phase_usage.phase_pos[Aqua]];
const double sw = std::max(s[Aqua], swco);
const double sg = s[Vapour];
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
// xw and xg are the fractions occupied by water and gas zones.
const double ssw = sg + sw;
// d = ssg = sw - swco + sg (using 'd' for consistency with mrst docs).
const double d = sg + sw - swco;
const double xw = (sw - swco) / d;
const double krw = krw_(sw);
const double krg = krg_(sg);
const double krow = krow_(ssw);
const double krog = krog_(d);
const double xg = 1 - xw;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = krw_.derivative(sw);
const double dkrgg = krg_.derivative(sg);
const double dkrow = krow_.derivative(ssw);
const double dkrog = krog_.derivative(d);
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Liquid + Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[Liquid + Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[Vapour + Vapour*np] = dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncGwsegNonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncGwsegNonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
} // namespace Opm
// SatFuncGwseg.cpp can now be removed
} // namespace Opm

643
opm/core/props/satfunc/SatFuncGwseg.hpp
View File

@@ -19,81 +19,602 @@
#ifndef SATFUNCGWSEG_HPP
#define SATFUNCGWSEG_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
class SatFuncGwsegUniform : public BlackoilPhases
template<class TableType>
class SatFuncGwseg : public SatFuncBase<TableType>
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
void evalKr(const double* s, double* kr, const EPSTransforms* epst) const;
void evalKr(const double* s, double* kr, 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;
void evalPc(const double* s, double* pc, const EPSTransforms* epst) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const;
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
typedef SatFuncGwseg<UniformTableLinear<double> > SatFuncGwsegUniform;
typedef SatFuncGwseg<NonuniformTableLinear<double> > SatFuncGwsegNonuniform;
class SatFuncGwsegNonuniform : public BlackoilPhases
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr) const
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
if (this->phase_usage.num_phases == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
const double sw = std::max(s[BlackoilPhases::Aqua], swco);
const double sg = s[BlackoilPhases::Vapour];
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
// xw and xg are the fractions occupied by water and gas zones.
const double ssg = sw - swco + sg;
const double xw = (sw - swco) / ssg;
const double xg = 1 - xw;
const double ssw = sg + sw;
const double krw = this->krw_(sw);
const double krg = this->krg_(sg);
const double krow = this->krow_(ssw);
const double krog = this->krog_(ssg);
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr, const EPSTransforms* epst) const
{
if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco);
const double sg = s[BlackoilPhases::Vapour];
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKr(const double* s, double* kr, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{
if (this->phase_usage.num_phases == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco);
const double sg = s[BlackoilPhases::Vapour];
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
// The code below corresponds to EHYSTR * 0 * * KR/
// - wettability properties water>oil>gas.
// - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase.
// The imb-curve currently only differs from drainage curves via endpoint scaling ...
// Water - use drainage curve only
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
// Gas
double krg;
if (sg >= sat_hyst->sg_hyst) { // Drainage
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
} else { // Imbibition
double sg_shifted = sg + sat_hyst->sg_shift;
double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_);
}
// Oil in water
double krow;
if (ssow >= sat_hyst->sow_hyst) { // Drainage
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
} else { // Imbibition
double ssow_shifted = ssow + sat_hyst->sow_shift;
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_);
}
// Oil in gas and connate water - use drainage curve only
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
// relperms
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
double swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
const double sw = std::max(s[BlackoilPhases::Aqua], swco);
const double sg = s[BlackoilPhases::Vapour];
const double eps = 1e-5;
swco = std::min(swco,sw-eps);
// xw and xg are the fractions occupied by water and gas zones.
const double ssw = sg + sw;
// d = ssg = sw - swco + sg (using 'd' for consistency with mrst docs).
const double d = sg + sw - swco;
const double xw = (sw - swco) / d;
const double krw = this->krw_(sw);
const double krg = this->krg_(sg);
const double krow = this->krow_(ssw);
const double krog = this->krog_(d);
const double xg = 1 - xw;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
// Derivatives.
const double dkrww = this->krw_.derivative(sw);
const double dkrgg = this->krg_.derivative(sg);
const double dkrow = this->krow_.derivative(ssw);
const double dkrog = this->krog_.derivative(d);
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco);
const double sg = s[BlackoilPhases::Vapour];
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
const double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
const double krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
const double krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
const double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
// Derivatives.
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
double dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
double dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst, const EPSTransforms* epst_hyst, const SatHyst* sat_hyst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// Relative permeability model based on segregation of water
// and gas, with oil present in both water and gas zones.
// TODO Also consider connate gas ...
double _swco = this->smin_[this->phase_usage.phase_pos[BlackoilPhases::Aqua]];
double swco = epst->wat.smin;
const double sw = std::max(s[BlackoilPhases::Aqua], swco);
const double sg = s[BlackoilPhases::Vapour];
const double eps = 1e-6;
swco = std::min(swco,sw-eps);
const double ssw = sg + sw;
const double ssg = std::max(sg + sw - swco, eps);
const double d = ssg; // = sw - swco + sg (using 'd' for consistency with mrst docs).
double ssow = 1.0-ssw;
double ssog = 1.0-ssg-swco;
// The code below corresponds to EHYSTR * 0 * * KR/
// - wettability properties water>oil>gas.
// - Carlsen hysteresis model for non-wetting (scanning=shifted_imb). No hysteresis for wetting phase.
// The imb-curve currently only differs from drainage curves via endpoint scaling ...
// Water - use drainage curve only
double _sw = epst->wat.scaleSat(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(sw, 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double krw = epst->wat.scaleKr(sw, this->krw_(_sw), this->krwr_);
double dkrww = _dsdsw*epst->wat.scaleKrDeriv(sw, this->krw_.derivative(_sw));
// Gas
double krg, dkrgg;
if (sg >= sat_hyst->sg_hyst) { // Drainage
double _sg = epst->gas.scaleSat(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(sg, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst->gas.scaleKr(sg, this->krg_(_sg), this->krgr_);
dkrgg = _dsdsg*epst->gas.scaleKrDeriv(sg, this->krg_.derivative(_sg));
} else { // Imbibition
double sg_shifted = sg + sat_hyst->sg_shift;
double _sg = epst_hyst->gas.scaleSat(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst_hyst->gas.scaleSatDeriv(sg_shifted, 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
krg = epst_hyst->gas.scaleKr(sg_shifted, this->krg_(_sg), this->krgr_);
dkrgg = _dsdsg*epst_hyst->gas.scaleKrDeriv(sg_shifted, this->krg_.derivative(_sg));
}
// Oil in water
double krow, dkrow;
if (ssow >= sat_hyst->sow_hyst) { // Drainage
double _ssow = epst->watoil.scaleSat(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst->watoil.scaleSatDeriv(ssow, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst->watoil.scaleKr(ssow, this->krow_(1.0-_ssow), this->krorw_);
dkrow = _dsdssow*epst->watoil.scaleKrDeriv(ssow, this->krow_.derivative(1.0-_ssow));
} else { // Imbibition
double ssow_shifted = ssow + sat_hyst->sow_shift;
double _ssow = epst_hyst->watoil.scaleSat(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssow = epst_hyst->watoil.scaleSatDeriv(ssow_shifted, 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
krow = epst_hyst->watoil.scaleKr(ssow_shifted, this->krow_(1.0-_ssow), this->krorw_);
dkrow = _dsdssow*epst_hyst->watoil.scaleKrDeriv(ssow_shifted, this->krow_.derivative(1.0-_ssow));
}
// Oil in gas and connate water - use drainage curve only
double _ssog = epst->gasoil.scaleSat(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdssog = epst->gasoil.scaleSatDeriv(ssog, 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krog = epst->gasoil.scaleKr(ssog, this->krog_(1.0-_ssog-_swco), this->krorg_);
double dkrog = _dsdssog*epst->gasoil.scaleKrDeriv(ssog, this->krog_.derivative(1.0-_ssog-_swco));
// xw and xg are the fractions occupied by water and gas zones.
const double xw = (sw - swco) / d;
const double xg = 1 - xw;
// relperms
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = xw*krow + xg*krog;
// Derivatives.
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = (xg/d)*krow + xw*dkrow - (xg/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = -(xw/d)*krow + xw*dkrow + (xw/d)*krog + xg*dkrog;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
return;
}
OPM_THROW(std::runtime_error, "SatFuncGwseg -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPc(const double* s, double* pc, const EPSTransforms* epst) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
double _sw = epst->wat.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcow_(_sw);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcog_(_sg);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
}
}
template<class TableType>
void SatFuncGwseg<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds, const EPSTransforms* epst) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
double _sw = epst->wat.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcow_(s[pos]);
double _dsdsw = epst->wat.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = _dsdsw*this->pcow_.derivative(_sw);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
double _sg = epst->gas.scaleSatPc(s[pos], this->smin_[pos], this->smax_[pos]);
pc[pos] = this->pcog_(_sg);
double _dsdsg = epst->gas.scaleSatDerivPc(s[pos], this->smin_[pos], this->smax_[pos]);
dpcds[np*pos + pos] = _dsdsg*this->pcog_.derivative(_sg);
}
}
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
} // namespace Opm
#endif // SATFUNCGWSEG_HPP

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

@@ -30,536 +30,194 @@
namespace Opm
{
void SatFuncSimpleUniform::init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples)
// SatFuncSimple.cpp can now be removed and the code below moved to SatFuncBase.cpp ...
template<>
void SatFuncBase<NonuniformTableLinear<double> >::initializeTableType(NonuniformTableLinear<double> & table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples)
{
phase_usage = phase_usg;
double swco = 0.0;
double swmax = 1.0;
if (phase_usage.phase_used[Aqua]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
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);
SatFuncSimpleUniform::ExtendTable(sw,sw_ex,1);
SatFuncSimpleUniform::ExtendTable(krw,krw_ex,0);
SatFuncSimpleUniform::ExtendTable(krow,krow_ex,0);
SatFuncSimpleUniform::ExtendTable(pcow,pcow_ex,0);
buildUniformMonotoneTable(sw_ex, krw_ex, samples, krw_);
buildUniformMonotoneTable(sw_ex, krow_ex, samples, krow_);
buildUniformMonotoneTable(sw_ex, pcow_ex, samples, pcow_);
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]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncSimpleUniform::ExtendTable(sg,sg_ex,1);
SatFuncSimpleUniform::ExtendTable(krg,krg_ex,0);
SatFuncSimpleUniform::ExtendTable(krog,krog_ex,0);
SatFuncSimpleUniform::ExtendTable(pcog,pcog_ex,0);
buildUniformMonotoneTable(sg_ex, krg_ex, samples, krg_);
buildUniformMonotoneTable(sg_ex, krog_ex, samples, krog_);
buildUniformMonotoneTable(sg_ex, pcog_ex, samples, pcog_);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
table = NonuniformTableLinear<double>(arg, value);
}
void SatFuncSimpleUniform::ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
template<>
void SatFuncBase<UniformTableLinear<double> >::initializeTableType(UniformTableLinear<double> & table,
const std::vector<double>& arg,
const std::vector<double>& value,
const int samples)
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
buildUniformMonotoneTable(arg, value, samples, table);
}
void SatFuncSimpleUniform::evalKr(const double* s, double* kr) const
double EPSTransforms::Transform::scaleSat(double s, double s_r, double s_cr, double s_max) const
{
if (phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
if (doNotScale) {
return s;
} else if (!do_3pt) { // 2-pt
if (s <= scr) {
return s_cr;
} else {
return (s >= smax) ? s_max : s_cr + (s-scr)*slope1;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double so = s[opos];
double krow = krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
} else if (s <= sr) {
return (s <= scr) ? s_cr : s_cr+(s-scr)*slope1;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
return (s >= smax) ? s_max : s_r+(s-sr)*slope2;
}
}
void SatFuncSimpleUniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
double EPSTransforms::Transform::scaleSatInv(double ss, double s_r, double s_cr, double s_max) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
if (doNotScale) {
return ss;
} else if (!do_3pt) { // 2-pt
if (ss <= s_cr) {
return scr;
} else {
return (ss >= s_max) ? smax : scr + (ss-s_cr)/slope1;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[Liquid + Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + Vapour*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double so = s[opos];
double krow = krow_(1.0-so);
double dkrow = krow_.derivative(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else if (ss <= s_r) {
return (ss <= s_cr) ? scr : scr+(ss-s_cr)/slope1;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncSimpleUniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
return (ss >= s_max) ? smax : sr+(ss-s_r)/slope2;
}
}
void SatFuncSimpleUniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
double EPSTransforms::Transform::scaleSatDeriv(double s, double s_r, double s_cr, double s_max) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncSimpleNonuniform ======
void SatFuncSimpleNonuniform::ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const
{
int n = xv.size();
xv_ex[0] = xv[0]-pm;
xv_ex[n+1] = xv[n-1]+pm;
for (int i=0; i<n; i++)
{
xv_ex[i+1] = xv[i];
}
}
void SatFuncSimpleNonuniform::init(const EclipseGridParser& deck,
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]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
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)");
if (doNotScale) {
return 1.0;
} else if (!do_3pt) { // 2-pt
if (s <= scr) {
return 0.0;
} else {
return (s >= smax) ? 0.0 : slope1;
}
// 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);
SatFuncSimpleNonuniform::ExtendTable(sw,sw_ex,1);
SatFuncSimpleNonuniform::ExtendTable(krw,krw_ex,0);
SatFuncSimpleNonuniform::ExtendTable(krow,krow_ex,0);
SatFuncSimpleNonuniform::ExtendTable(pcow,pcow_ex,0);
krw_ = NonuniformTableLinear<double>(sw_ex, krw_ex);
krow_ = NonuniformTableLinear<double>(sw_ex, krow_ex);
pcow_ = NonuniformTableLinear<double>(sw_ex, pcow_ex);
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]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
// 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);
SatFuncSimpleNonuniform::ExtendTable(sg,sg_ex,1);
SatFuncSimpleNonuniform::ExtendTable(krg,krg_ex,0);
SatFuncSimpleNonuniform::ExtendTable(krog,krog_ex,0);
SatFuncSimpleNonuniform::ExtendTable(pcog,pcog_ex,0);
krg_ = NonuniformTableLinear<double>(sg_ex, krg_ex);
krog_ = NonuniformTableLinear<double>(sg_ex, krog_ex);
pcog_ = NonuniformTableLinear<double>(sg_ex, pcog_ex);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncSimpleNonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double so = s[opos];
double krow = krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
} else if (s <= sr) {
return (s <= scr) ? 0.0 : slope1;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
return (s >= smax) ? 0.0 : slope2;
}
}
void SatFuncSimpleNonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
double EPSTransforms::Transform::scaleSatPc(double s, double s_min, double s_max) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[Liquid + Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + Vapour*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double so = s[opos];
double krow = krow_(1.0-so);
double dkrow = krow_.derivative(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
if (doNotScale) {
return s;
} else if (s<=smin) {
return s_min;
} else if (s <= smax) {
return s_min + (s-smin)*(s_max-s_min)/(smax-smin);
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
return s_max;
}
}
void SatFuncSimpleNonuniform::evalPc(const double* s, double* pc) const
double EPSTransforms::Transform::scaleSatDerivPc(double s, double s_min, double s_max) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
if (doNotScale) {
return 1.0;
} else if (s<smin) {
return 0.0;
} else if (s <= smax) {
return (s_max-s_min)/(smax-smin);
} else {
return 0.0;
}
}
void SatFuncSimpleNonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
double EPSTransforms::Transform::scaleKr(double s, double kr, double krsr_tab) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
if (doKrCrit) {
if (s <= scr) {
return 0.0;
} else if (s <= sr) {
return kr*krSlopeCrit;
} else if (s <= smax) {
if (doSatInterp)
return krsr + (s-sr)*krSlopeMax; // Note: Scaling independent of kr-value ...
else
return krsr + (kr-krsr_tab)*krSlopeMax;
} else {
return krmax;
}
} else if (doKrMax) {
if (s <= scr) {
return 0.0;
} else if (s <= smax) {
return kr*krSlopeMax;
} else {
return krmax;
}
} else {
return kr;
}
}
double EPSTransforms::Transform::scaleKrDeriv(double s, double krDeriv) const
{
if (doKrCrit) {
if (s <= scr) {
return 0.0;
} else if (s <= sr) {
return krDeriv*krSlopeCrit;
} else if (s <= smax) {
if (doSatInterp)
return krSlopeMax; // Note: Scaling independent of kr-value ...
else
return krDeriv*krSlopeMax;
} else {
return 0.0;
}
} else if (doKrMax) {
if (s <= scr) {
return 0.0;
} else if (s <= smax) {
return krDeriv*krSlopeMax;
} else {
return 0.0;
}
} else {
if (s <= scr) {
return 0.0;
} else if (s <= smax) {
return krDeriv;
} else {
return 0.0;
}
}
}
void EPSTransforms::Transform::printMe(std::ostream & out)
{
out << "doNotScale: " << doNotScale << std::endl;
out << "do_3pt: " << do_3pt << std::endl;
out << "smin: " << smin << std::endl;
out << "scr: " << scr << std::endl;
out << "sr: " << sr << std::endl;
out << "smax: " << smax << std::endl;
out << "slope1: " << slope1 << std::endl;
out << "slope2: " << slope2 << std::endl;
out << "doKrMax: " << doKrMax << std::endl;
out << "doKrCrit: " << doKrCrit << std::endl;
out << "doSatInterp: " << doSatInterp << std::endl;
out << "krsr: " << krsr << std::endl;
out << "krmax: " << krmax << std::endl;
out << "krSlopeMax: " << krSlopeMax << std::endl;
out << "krSlopeCrit: " << krSlopeCrit << std::endl;
}
void SatHyst::printMe(std::ostream & out)
{
out << "sg_hyst: " << sg_hyst << std::endl;
out << "sg_shift: " << sg_shift << std::endl;
out << "sow_hyst: " << sow_hyst << std::endl;
out << "sow_shift: " << sow_shift << std::endl;
};
} // namespace Opm

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

@@ -19,80 +19,264 @@
#ifndef SATFUNCSIMPLE_HPP
#define SATFUNCSIMPLE_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
class SatFuncSimpleUniform : public BlackoilPhases
template<class TableType>
class SatFuncSimple : public SatFuncBase<TableType>
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
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
{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
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
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
{OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");}
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
typedef SatFuncSimple<UniformTableLinear<double> > SatFuncSimpleUniform;
typedef SatFuncSimple<NonuniformTableLinear<double> > SatFuncSimpleNonuniform;
class SatFuncSimpleNonuniform : public BlackoilPhases
template<class TableType>
void SatFuncSimple<TableType>::evalKr(const double* s, double* kr) const
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
void ExtendTable(const std::vector<double>& xv,
std::vector<double>& xv_ex,
double pm) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
if (this->phase_usage.num_phases == 3) {
// A simplified relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double krw = this->krw_(sw);
double krg = this->krg_(sg);
double krow = this->krow_(sw + sg); // = 1 - so
// double krog = krog_(sg); // = 1 - so - sw
// double krocw = krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krow;
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double so = s[opos];
double krow = this->krow_(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// A simplified relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krow = this->krow_(sw + sg);
double dkrow = this->krow_.derivative(sw + sg);
// double krog = krog_(sg);
// double dkrog = krog_.derivative(sg);
// double krocw = krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = s[opos];
double krow = this->krow_(1.0-so);
double dkrow = this->krow_.derivative(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds, const EPSTransforms* epst) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
// A simplified relative permeability model.
// Define KR(s) = scaleKr(kr(scalSat(s)))
// Thus KR'(s) = scaleKr'(kr(scaleSat(s)))*kr'((scaleSat(s))*scaleSat'(s)
double _sw = epst->wat.scaleSat(s[wpos], 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _dsdsw = epst->wat.scaleSatDeriv(s[wpos], 1.0-this->sowcr_-this->smin_[gpos], this->swcr_, this->smax_[wpos]);
double _sg = epst->gas.scaleSat(s[gpos], 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _dsdsg = epst->gas.scaleSatDeriv(s[gpos], 1.0-this->sogcr_-this->smin_[wpos], this->sgcr_, this->smax_[gpos]);
double _sow = epst->watoil.scaleSat(1.0-s[wpos]-s[gpos], 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double _dsdsow = epst->watoil.scaleSatDeriv(1.0-s[wpos]-s[gpos], 1.0-this->swcr_-this->smin_[gpos], this->sowcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
//double _sog = epst->gasoil.scaleSat(1.0-s[wpos]-s[gpos], 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
//double _dsdsog = epst->gasoil.scaleSatDeriv(1.0-s[wpos]-s[gpos], 1.0-this->sgcr_-this->smin_[wpos], this->sogcr_, 1.0-this->smin_[wpos]-this->smin_[gpos]);
double krw = epst->wat.scaleKr(s[wpos], this->krw_(_sw), this->krwr_);
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...
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_); // ????
//double dkrog = _dsdsog*epst->gasoil.scaleKrDeriv(1.0-s[wpos]-s[gpos], this->krog_.derivative(1.0-_sog-this->smin_[wpos])); // ????
// double krocw = krocw_;
kr[wpos] = krw;
kr[gpos] = krg;
kr[BlackoilPhases::Liquid] = krow;
//krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
dkrds[wpos + wpos*np] = dkrww;
dkrds[gpos + gpos*np] = dkrgg;
//dkrds[Liquid + Aqua*np] = dkrow;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Liquid*np] = -dkrow;
//krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[BlackoilPhases::Liquid + gpos*np] = 0.0;
//krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
//+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
OPM_THROW(std::runtime_error, "SatFuncSimple -- need to be implemented ...");
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double so = s[opos];
double krow = this->krow_(1.0-so);
double dkrow = this->krow_.derivative(1.0-so);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalPc(const double* s, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
}
}
template<class TableType>
void SatFuncSimple<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
}
}
} // namespace Opm
#endif // SATFUNCSIMPLE_HPP

389
opm/core/props/satfunc/SatFuncStone2.cpp
View File

@@ -29,390 +29,5 @@
namespace Opm
{
void SatFuncStone2Uniform::init(const EclipseGridParser& deck,
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]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
buildUniformMonotoneTable(sw, krw, samples, krw_);
buildUniformMonotoneTable(sw, krow, samples, krow_);
buildUniformMonotoneTable(sw, pcow, samples, pcow_);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
buildUniformMonotoneTable(sg, krg, samples, krg_);
buildUniformMonotoneTable(sg, krog, samples, krog_);
buildUniformMonotoneTable(sg, pcog, samples, pcog_);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncStone2Uniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
double krog = krog_(sg); // = 1 - so - sw
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncStone2Uniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncStone2Uniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncStone2Uniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
// ====== Methods for SatFuncStone2Nonuniform ======
void SatFuncStone2Nonuniform::init(const EclipseGridParser& deck,
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]) {
const SWOF::table_t& swof_table = deck.getSWOF().swof_;
const std::vector<double>& sw = swof_table[table_num][0];
const std::vector<double>& krw = swof_table[table_num][1];
const std::vector<double>& krow = swof_table[table_num][2];
const std::vector<double>& pcow = swof_table[table_num][3];
krw_ = NonuniformTableLinear<double>(sw, krw);
krow_ = NonuniformTableLinear<double>(sw, krow);
pcow_ = NonuniformTableLinear<double>(sw, pcow);
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();
}
if (phase_usage.phase_used[Vapour]) {
const SGOF::table_t& sgof_table = deck.getSGOF().sgof_;
const std::vector<double>& sg = sgof_table[table_num][0];
const std::vector<double>& krg = sgof_table[table_num][1];
const std::vector<double>& krog = sgof_table[table_num][2];
const std::vector<double>& pcog = sgof_table[table_num][3];
krg_ = NonuniformTableLinear<double>(sg, krg);
krog_ = NonuniformTableLinear<double>(sg, krog);
pcog_ = NonuniformTableLinear<double>(sg, pcog);
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();
}
// These only consider water min/max sats. Consider gas sats?
smin_[phase_usage.phase_pos[Liquid]] = 1.0 - swmax;
smax_[phase_usage.phase_pos[Liquid]] = 1.0 - swco;
}
void SatFuncStone2Nonuniform::evalKr(const double* s, double* kr) const
{
if (phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double krg = krg_(sg);
double krow = krow_(sw + sg); // = 1 - so
double krog = krog_(sg); // = 1 - so - sw
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double krow = krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double krog = krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
void SatFuncStone2Nonuniform::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[Aqua];
double sg = s[Vapour];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krow = krow_(sw + sg);
double dkrow = krow_.derivative(sw + sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
double krocw = krocw_;
kr[Aqua] = krw;
kr[Vapour] = krg;
kr[Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[Liquid] < 0.0) {
kr[Liquid] = 0.0;
}
dkrds[Aqua + Aqua*np] = dkrww;
dkrds[Vapour + Vapour*np] = dkrgg;
dkrds[Liquid + Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[Liquid + Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (phase_usage.phase_used[Aqua]) {
int wpos = phase_usage.phase_pos[Aqua];
int opos = phase_usage.phase_pos[Liquid];
double sw = s[wpos];
double krw = krw_(sw);
double dkrww = krw_.derivative(sw);
double krow = krow_(sw);
double dkrow = krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(phase_usage.phase_used[Vapour]);
int gpos = phase_usage.phase_pos[Vapour];
int opos = phase_usage.phase_pos[Liquid];
double sg = s[gpos];
double krg = krg_(sg);
double dkrgg = krg_.derivative(sg);
double krog = krog_(sg);
double dkrog = krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
void SatFuncStone2Nonuniform::evalPc(const double* s, double* pc) const
{
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
}
}
void SatFuncStone2Nonuniform::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[phase_usage.phase_pos[Liquid]] = 0.0;
if (phase_usage.phase_used[Aqua]) {
int pos = phase_usage.phase_pos[Aqua];
pc[pos] = pcow_(s[pos]);
dpcds[np*pos + pos] = pcow_.derivative(s[pos]);
}
if (phase_usage.phase_used[Vapour]) {
int pos = phase_usage.phase_pos[Vapour];
pc[pos] = pcog_(s[pos]);
dpcds[np*pos + pos] = pcog_.derivative(s[pos]);
}
}
} // namespace Opm
// SatFuncStone2.cpp can now be removed
} // namespace Opm

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

@@ -19,74 +19,180 @@
#ifndef SATFUNCSTONE2_HPP
#define SATFUNCSTONE2_HPP
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/UniformTableLinear.hpp>
#include <opm/core/utility/NonuniformTableLinear.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <vector>
#include <opm/core/props/satfunc/SatFuncBase.hpp>
namespace Opm
{
class SatFuncStone2Uniform : public BlackoilPhases
template<class TableType>
class SatFuncStone2 : public SatFuncBase<TableType>
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
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
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
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
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
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
{OPM_THROW(std::runtime_error, "SatFuncStone2 -- need to be implemented ...");}
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
UniformTableLinear<double> krw_;
UniformTableLinear<double> krow_;
UniformTableLinear<double> pcow_;
UniformTableLinear<double> krg_;
UniformTableLinear<double> krog_;
UniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
typedef SatFuncStone2<UniformTableLinear<double> > SatFuncStone2Uniform;
typedef SatFuncStone2<NonuniformTableLinear<double> > SatFuncStone2Nonuniform;
class SatFuncStone2Nonuniform : public BlackoilPhases
template<class TableType>
void SatFuncStone2<TableType>::evalKr(const double* s, double* kr) const
{
public:
void init(const EclipseGridParser& deck,
const int table_num,
const PhaseUsage phase_usg,
const int samples);
void evalKr(const double* s, double* kr) const;
void evalKrDeriv(const double* s, double* kr, double* dkrds) const;
void evalPc(const double* s, double* pc) const;
void evalPcDeriv(const double* s, double* pc, double* dpcds) const;
double smin_[PhaseUsage::MaxNumPhases];
double smax_[PhaseUsage::MaxNumPhases];
double krwmax_; // Max water relperm
double kromax_; // Max oil relperm
double swcr_; // Critical water saturation.
double krwr_; // Water relperm at critical oil-in-water saturation.
double sowcr_; // Critical oil-in-water saturation.
double krorw_; // Oil relperm at critical water saturation.
private:
PhaseUsage phase_usage; // A copy of the outer class' phase_usage_.
NonuniformTableLinear<double> krw_;
NonuniformTableLinear<double> krow_;
NonuniformTableLinear<double> pcow_;
NonuniformTableLinear<double> krg_;
NonuniformTableLinear<double> krog_;
NonuniformTableLinear<double> pcog_;
double krocw_; // = krow_(s_wc)
};
if (this->phase_usage.num_phases == 3) {
// Stone-II relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double krw = this->krw_(sw);
double krg = this->krg_(sg);
double krow = this->krow_(sw + sg); // = 1 - so
double krog = this->krog_(sg); // = 1 - so - sw
double krocw = this->krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double krow = this->krow_(sw);
kr[wpos] = krw;
kr[opos] = krow;
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double krog = this->krog_(sg);
kr[gpos] = krg;
kr[opos] = krog;
}
}
template<class TableType>
void SatFuncStone2<TableType>::evalKrDeriv(const double* s, double* kr, double* dkrds) const
{
const int np = this->phase_usage.num_phases;
std::fill(dkrds, dkrds + np*np, 0.0);
if (np == 3) {
// Stone-II relative permeability model.
double sw = s[BlackoilPhases::Aqua];
double sg = s[BlackoilPhases::Vapour];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krow = this->krow_(sw + sg);
double dkrow = this->krow_.derivative(sw + sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
double krocw = this->krocw_;
kr[BlackoilPhases::Aqua] = krw;
kr[BlackoilPhases::Vapour] = krg;
kr[BlackoilPhases::Liquid] = krocw*((krow/krocw + krw)*(krog/krocw + krg) - krw - krg);
if (kr[BlackoilPhases::Liquid] < 0.0) {
kr[BlackoilPhases::Liquid] = 0.0;
}
dkrds[BlackoilPhases::Aqua + BlackoilPhases::Aqua*np] = dkrww;
dkrds[BlackoilPhases::Vapour + BlackoilPhases::Vapour*np] = dkrgg;
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Aqua*np] = krocw*((dkrow/krocw + dkrww)*(krog/krocw + krg) - dkrww);
dkrds[BlackoilPhases::Liquid + BlackoilPhases::Vapour*np] = krocw*((krow/krocw + krw)*(dkrog/krocw + dkrgg) - dkrgg)
+ krocw*((dkrow/krocw + krw)*(krog/krocw + krg) - dkrgg);
return;
}
// We have a two-phase situation. We know that oil is active.
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int wpos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sw = s[wpos];
double krw = this->krw_(sw);
double dkrww = this->krw_.derivative(sw);
double krow = this->krow_(sw);
double dkrow = this->krow_.derivative(sw);
kr[wpos] = krw;
kr[opos] = krow;
dkrds[wpos + wpos*np] = dkrww;
dkrds[opos + wpos*np] = dkrow; // Row opos, column wpos, fortran order.
} else {
assert(this->phase_usage.phase_used[BlackoilPhases::Vapour]);
int gpos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
int opos = this->phase_usage.phase_pos[BlackoilPhases::Liquid];
double sg = s[gpos];
double krg = this->krg_(sg);
double dkrgg = this->krg_.derivative(sg);
double krog = this->krog_(sg);
double dkrog = this->krog_.derivative(sg);
kr[gpos] = krg;
kr[opos] = krog;
dkrds[gpos + gpos*np] = dkrgg;
dkrds[opos + gpos*np] = dkrog;
}
}
template<class TableType>
void SatFuncStone2<TableType>::evalPc(const double* s, double* pc) const
{
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
}
}
template<class TableType>
void SatFuncStone2<TableType>::evalPcDeriv(const double* s, double* pc, double* dpcds) const
{
// The problem of determining three-phase capillary pressures
// is very hard experimentally, usually one extends two-phase
// data (as for relative permeability).
// In our approach the derivative matrix is quite sparse, only
// the diagonal elements corresponding to non-oil phases are
// (potentially) nonzero.
const int np = this->phase_usage.num_phases;
std::fill(dpcds, dpcds + np*np, 0.0);
pc[this->phase_usage.phase_pos[BlackoilPhases::Liquid]] = 0.0;
if (this->phase_usage.phase_used[BlackoilPhases::Aqua]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Aqua];
pc[pos] = this->pcow_(s[pos]);
dpcds[np*pos + pos] = this->pcow_.derivative(s[pos]);
}
if (this->phase_usage.phase_used[BlackoilPhases::Vapour]) {
int pos = this->phase_usage.phase_pos[BlackoilPhases::Vapour];
pc[pos] = this->pcog_(s[pos]);
dpcds[np*pos + pos] = this->pcog_.derivative(s[pos]);
}
}
} // namespace Opm
#endif // SATFUNCSTONE2_HPP

118
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;
@@ -82,6 +85,38 @@ namespace Opm
int dimensions,
const int samples);
/// Initialize from deck and grid.
/// \param[in] newParserDeck 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);
/// Initialize from deck and grid.
/// \param[in] newParserDeck Deck input parser
/// \param[in] number_of_cells The number of cells of the 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] global_cell The mapping from local cell indices of the grid to
/// global cell indices used in the deck.
/// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid.
/// \param[in] dimensions The dimensions of the grid.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples);
/// \return P, the number of phases.
int numPhases() const;
@@ -123,37 +158,92 @@ namespace Opm
const int* cells,
double* smin,
double* smax) const;
/// Update saturation state for the hysteresis tracking
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
void updateSatHyst(const int n,
const int* cells,
const double* s);
private:
PhaseUsage phase_usage_;
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
struct { // End point scaling parameters
std::vector<double> swl_;
std::vector<double> swcr_;
std::vector<double> swu_;
std::vector<double> sowcr_;
std::vector<double> krw_;
std::vector<double> krwr_;
std::vector<double> kro_;
std::vector<double> krorw_;
} eps_;
bool do_eps_; // ENDSCALE is active
bool do_3pt_; // SCALECRS: YES~true NO~false
bool do_hyst_; // Keywords ISWL etc detected
std::vector<EPSTransforms> eps_transf_;
std::vector<EPSTransforms> eps_transf_hyst_;
std::vector<SatHyst> sat_hyst_;
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const;
int dimensions);
template<class T>
void initEPSHyst(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
template<class T>
void initEPS(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
void initEPSParam(const int cell,
EPSTransforms::Transform& data,
const bool oil,
const double sl_tab,
const double scr_tab,
const double su_tab,
const double sxcr_tab,
const double s0_tab,
const double krsr_tab,
const double krmax_tab,
const std::vector<double>& sl,
const std::vector<double>& scr,
const std::vector<double>& su,
const std::vector<double>& sxcr,
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; }
};

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

File diff suppressed because it is too large Load Diff

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

@@ -33,6 +33,8 @@
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
#include <iostream>
#include <cmath>
@@ -547,6 +549,103 @@ 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")) {
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,
@@ -666,10 +765,6 @@ namespace Opm
begin_cell_centroids, 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.
@@ -904,6 +999,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

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

@@ -68,16 +68,18 @@ namespace Opm
{
static double handleBracketingFailure(const double x0, const double x1, const double f0, const double f1)
{
OPM_MESSAGE("Error in parameters, zero not bracketed: [a, b] = ["
<< x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1
<< "");
OPM_REPORT;
std::cerr << "Error in parameters, zero not bracketed: [a, b] = ["
<< x0 << ", " << x1 << "] f(a) = " << f0 << " f(b) = " << f1
<< "";
return std::fabs(f0) < std::fabs(f1) ? x0 : x1;
}
static double handleTooManyIterations(const double x0, const double x1, const int maxiter)
{
OPM_MESSAGE("Maximum number of iterations exceeded: " << maxiter
<< ", current interval is [" << std::min(x0, x1) << ", "
<< std::max(x0, x1) << "]");
OPM_REPORT;
std::cerr << "Maximum number of iterations exceeded: " << maxiter
<< ", current interval is [" << std::min(x0, x1) << ", "
<< std::max(x0, x1) << "]";
return 0.5*(x0 + x1);
}
};

9
opm/core/utility/Spline.hpp
View File

@@ -945,6 +945,7 @@ public:
int r = 3;
if (x0 < xMin()) {
static_cast<void>(extrapolate);
assert(extrapolate);
Scalar m = evalDerivative_(xMin(), /*segmentIdx=*/0);
if (std::abs(m) < 1e-20)
@@ -1575,16 +1576,16 @@ protected:
{ return 2*3*t - 2; }
// third derivative of the hermite basis functions
Scalar h00_prime3_(Scalar t) const
Scalar h00_prime3_(Scalar /*t*/) const
{ return 2*3*2; }
Scalar h10_prime3_(Scalar t) const
Scalar h10_prime3_(Scalar /*t*/) const
{ return 2*3; }
Scalar h01_prime3_(Scalar t) const
Scalar h01_prime3_(Scalar /*t*/) const
{ return -2*3*2; }
Scalar h11_prime3_(Scalar t) const
Scalar h11_prime3_(Scalar /*t*/) const
{ return 2*3; }
// returns the monotonicality of an interval of a spline segment

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.

52
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_;
}
@@ -1141,4 +1141,54 @@ 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)) {
injection_specification.BHP_limit_ = well->getBHPLimit(timeStep);
injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(well->getInjectorType(timeStep)));
injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(well->getInjectorControlMode(timeStep)));
injection_specification.surface_flow_max_rate_ = well->getSurfaceInjectionRate(timeStep);
injection_specification.reservoir_flow_max_rate_ = well->getReservoirInjectionRate(timeStep);
production_specification.guide_rate_ = 0.0; // We know we're not a producer
}
else if (well->isProducer(timeStep)) {
production_specification.BHP_limit_ = well->getBHPLimit(timeStep);
production_specification.reservoir_flow_max_rate_ = well->getResVRate(timeStep);
production_specification.oil_max_rate_ = well->getOilRate(timeStep);
production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(well->getProducerControlMode(timeStep)));
production_specification.water_max_rate_ = well->getWaterRate(timeStep);
injection_specification.guide_rate_ = 0.0; // we know we're not an injector
}
std::shared_ptr<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 */

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

@@ -277,7 +277,6 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
/// Construct wells from deck.
WellsManager::WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability,
bool checkCellExistence)
@@ -314,86 +313,20 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
well_names.reserve(wells.size());
well_data.reserve(wells.size());
createWellsFromSpecs(wells, timeStep, grid, well_names, well_data, well_names_to_index, pu, cartesian_to_compressed, permeability);
setupWellControls(wells, timeStep, well_names, pu);
if (deck.hasField("WELOPEN")) {
const WELOPEN& welopen = deck.getWELOPEN();
for (size_t i = 0; i < welopen.welopen.size(); ++i) {
WelopenLine line = welopen.welopen[i];
std::string wellname = line.well_;
std::map<std::string, int>::const_iterator it = well_names_to_index.find(wellname);
if (it == well_names_to_index.end()) {
OPM_THROW(std::runtime_error, "Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
}
const int index = it->second;
if (line.openshutflag_ == "SHUT") {
int cur_ctrl = well_controls_get_current(w_->ctrls[index]);
if (cur_ctrl >= 0) {
well_controls_invert_current(w_->ctrls[index]);
}
assert(well_controls_get_current(w_->ctrls[index]) < 0);
} else if (line.openshutflag_ == "OPEN") {
int cur_ctrl = well_controls_get_current(w_->ctrls[index]);
if (cur_ctrl < 0) {
well_controls_invert_current(w_->ctrls[index]);
}
assert(well_controls_get_current(w_->ctrls[index]) >= 0);
} else {
OPM_THROW(std::runtime_error, "Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
}
}
{
GroupTreeNodeConstPtr fieldNode = eclipseState->getSchedule()->getGroupTree(timeStep)->getNode("FIELD");
GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(fieldNode->name());
well_collection_.addField(fieldGroup, timeStep, pu);
addChildGroups(fieldNode, eclipseState->getSchedule(), timeStep, pu);
}
// Build the well_collection_ well group hierarchy.
if (deck.hasField("GRUPTREE")) {
std::cout << "Found gruptree" << std::endl;
const GRUPTREE& gruptree = deck.getGRUPTREE();
std::map<std::string, std::string>::const_iterator it = gruptree.tree.begin();
for( ; it != gruptree.tree.end(); ++it) {
well_collection_.addChild(it->first, it->second, deck);
}
}
for (auto wellIter = wells.begin(); wellIter != wells.end(); ++wellIter ) {
well_collection_.addChild((*wellIter)->name(), (*wellIter)->getGroupName(timeStep), deck);
well_collection_.addWell((*wellIter), timeStep, pu);
}
// Set the guide rates:
if (deck.hasField("WGRUPCON")) {
std::cout << "Found Wgrupcon" << std::endl;
WGRUPCON wgrupcon = deck.getWGRUPCON();
const std::vector<WgrupconLine>& lines = wgrupcon.wgrupcon;
std::cout << well_collection_.getLeafNodes().size() << std::endl;
for (size_t i = 0; i < lines.size(); i++) {
std::string name = lines[i].well_;
const int wix = well_names_to_index[name];
WellNode& wellnode = *well_collection_.getLeafNodes()[wix];
assert(wellnode.name() == name);
if (well_data[wix].type == PRODUCER) {
wellnode.prodSpec().guide_rate_ = lines[i].guide_rate_;
if (lines[i].phase_ == "OIL") {
wellnode.prodSpec().guide_rate_type_ = ProductionSpecification::OIL;
} else {
OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for producer "
<< name << " in WGRUPCON, cannot handle.");
}
} else if (well_data[wix].type == INJECTOR) {
wellnode.injSpec().guide_rate_ = lines[i].guide_rate_;
if (lines[i].phase_ == "RAT") {
wellnode.injSpec().guide_rate_type_ = InjectionSpecification::RAT;
} else {
OPM_THROW(std::runtime_error, "Guide rate type " << lines[i].phase_ << " specified for injector "
<< name << " in WGRUPCON, cannot handle.");
}
} else {
OPM_THROW(std::runtime_error, "Unknown well type " << well_data[wix].type << " for well " << name);
}
}
}
well_collection_.setWellsPointer(w_);
well_collection_.applyGroupControls();
@@ -1216,6 +1149,10 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
for (auto wellIter= wells.begin(); wellIter != wells.end(); ++wellIter) {
WellConstPtr well = (*wellIter);
if ( !( well->getStatus( timeStep ) == WellCommon::SHUT || well->getStatus( timeStep ) == WellCommon::OPEN) ) {
OPM_THROW(std::runtime_error, "Currently we do not support well status " << WellCommon::Status2String(well->getStatus( timeStep )));
}
if (well->isInjector(timeStep)) {
clear_well_controls(well_index, w_);
int ok = 1;
@@ -1427,4 +1364,13 @@ WellsManager::WellsManager(struct Wells* W, bool checkCellExistence)
}
}
void WellsManager::addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage) {
for (auto childIter = parentNode->begin(); childIter != parentNode->end(); ++childIter) {
GroupTreeNodeConstPtr childNode = (*childIter).second;
well_collection_.addGroup(schedule->getGroup(childNode->name()), parentNode->name(), timeStep, phaseUsage);
addChildGroups(childNode, schedule, timeStep, phaseUsage);
}
}
} // namespace Opm

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

@@ -25,6 +25,7 @@
#include <opm/core/wells/WellCollection.hpp>
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTree.hpp>
struct Wells;
struct UnstructuredGrid;
@@ -79,7 +80,6 @@ namespace Opm
WellsManager(const Opm::EclipseStateConstPtr eclipseState,
const size_t timeStep,
const Opm::EclipseGridParser& deck,
const UnstructuredGrid& grid,
const double* permeability,
bool checkCellExistence=true);
@@ -152,6 +152,8 @@ namespace Opm
const std::map<int,int> cartesian_to_compressed,
const double* permeability);
void addChildGroups(GroupTreeNodeConstPtr parentNode, ScheduleConstPtr schedule, size_t timeStep, const PhaseUsage& phaseUsage);
// Data

3
tests/testBlackoilState1.DATA
View File

@@ -13,3 +13,6 @@ DZV
ACTNUM
1 998*2 3 /
DEPTHZ
121*2000 /

3
tests/testBlackoilState2.DATA
View File

@@ -9,3 +9,6 @@ DYV
DZV
10*10 /
DEPTHZ
110*2000 /

2
tests/test_parser.cpp
View File

@@ -43,7 +43,7 @@ BOOST_AUTO_TEST_CASE(CreateParser)
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile( filename1 );
BOOST_CHECK_EQUAL( 5U , deck->size() );
BOOST_CHECK_EQUAL( 6U , deck->size() );
Opm::DeckItemConstPtr actnum = deck->getKeyword("ACTNUM")->getRecord(0)->getItem(0);
const std::vector<int>& actnum_data = actnum->getIntData();

22
tests/test_spline.cpp
View File

@@ -80,14 +80,14 @@ void testCommon(const Spline &sp,
// make sure the derivatives are consistent with the curve
int np = 3*n;
for (int i = 0; i < np; ++i) {
double x = sp.xMin() + (sp.xMax() - sp.xMin())*i/np;
double xval = sp.xMin() + (sp.xMax() - sp.xMin())*i/np;
// first derivative
double y1 = sp.eval(x+epsFD);
double y0 = sp.eval(x);
double y1 = sp.eval(xval+epsFD);
double y0 = sp.eval(xval);
double mFD = (y1 - y0)/epsFD;
double m = sp.evalDerivative(x);
double m = sp.evalDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error,
@@ -95,11 +95,11 @@ void testCommon(const Spline &sp,
" (" << mFD << " - " << m << " = " << mFD - m << ")!");
// second derivative
y1 = sp.evalDerivative(x+epsFD);
y0 = sp.evalDerivative(x);
y1 = sp.evalDerivative(xval+epsFD);
y0 = sp.evalDerivative(xval);
mFD = (y1 - y0)/epsFD;
m = sp.evalSecondDerivative(x);
m = sp.evalSecondDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error,
@@ -107,11 +107,11 @@ void testCommon(const Spline &sp,
" (" << mFD << " - " << m << " = " << mFD - m << ")!");
// Third derivative
y1 = sp.evalSecondDerivative(x+epsFD);
y0 = sp.evalSecondDerivative(x);
y1 = sp.evalSecondDerivative(xval+epsFD);
y0 = sp.evalSecondDerivative(xval);
mFD = (y1 - y0)/epsFD;
m = sp.evalThirdDerivative(x);
m = sp.evalThirdDerivative(xval);
if (std::abs( mFD - m ) > 1000*epsFD)
OPM_THROW(std::runtime_error,
@@ -331,7 +331,7 @@ void plot()
std::cout << "\n";
}
int main(int argc, char** argv)
int main()
{
try {
testAll();

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());
}

108
tests/test_wellsgroup.cpp Normal file
View File

@@ -0,0 +1,108 @@
/*
Copyright 2014 Statoil.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <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)) {
BOOST_CHECK_EQUAL(well->getSurfaceInjectionRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->injSpec().BHP_limit_);
BOOST_CHECK_EQUAL(well->getReservoirInjectionRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_);
}
if (well->isProducer(2)) {
BOOST_CHECK_EQUAL(well->getResVRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->prodSpec().BHP_limit_);
BOOST_CHECK_EQUAL(well->getOilRate(2), wellsGroup->prodSpec().oil_max_rate_);
BOOST_CHECK_EQUAL(well->getWaterRate(2), wellsGroup->prodSpec().water_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_);
}
}
}
BOOST_AUTO_TEST_CASE(ConstructGroupFromGroup) {
ParserPtr parser(new Parser());
boost::filesystem::path scheduleFile("wells_group.data");
DeckConstPtr deck = parser->parseFile(scheduleFile.string());
EclipseStateConstPtr eclipseState(new EclipseState(deck));
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
std::vector<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_);
}
}
}

31
tests/test_wellsmanager.cpp
View File

@@ -207,7 +207,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(eclipseState, 0, Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
std::cout << "Checking new well structure, epoch 0" << std::endl;
@@ -223,7 +223,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(eclipseState, 1,Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
std::cout << "Checking new well structure, epoch 1" << std::endl;
@@ -249,7 +249,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(eclipseState, 0, Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells(),false));
@@ -257,7 +257,15 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(eclipseState, 1,Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
}
Deck.setCurrentEpoch(2);
{
Opm::WellsManager wellsManager(eclipseState, 2, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
@@ -303,3 +311,18 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
}
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

@@ -16,6 +16,10 @@ DYV
DZV
10.0 20.0 30.0 10.0 5.0 /
DEPTHZ
121*2000 /
SCHEDULE
WELSPECS
@@ -30,11 +34,11 @@ COMPDAT
WCONPROD
'PROD1' 'OPEN' 'ORAT' 20000 4* 1000 /
/
/
WCONINJE
'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 /
/
/
DATES
@@ -49,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

18
tests/wells_manager_data_expanded.data
View File

@@ -16,6 +16,12 @@ DYV
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
@@ -66,7 +72,15 @@ WCONINJE
TSTEP
14.0
/
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