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

Merge remote-tracking branch 'upstream/master' into master-refactor-for-cpgrid-support

Browse Source 
Manually resolved conflicts in:
	opm/core/io/eclipse/EclipseWriter.cpp
	opm/core/io/eclipse/EclipseWriter.hpp
	opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
This commit is contained in:
Markus Blatt 2014-04-08 21:50:00 +02:00
commit ebc86bc624
36 changed files with 4847 additions and 689 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
CMakeLists.txt
View File

@ -47,9 +47,6 @@ include (CMakeLists_files.cmake)
macro (config_hook)
# opm_need_version_of ("dune-common")
opm_need_version_of ("dune-istl")
list (APPEND ${project}_CONFIG_IMPL_VARS
HAVE_DUNE_ISTL
)
endmacro (config_hook)
macro (prereqs_hook)
@ -74,6 +71,12 @@ macro (sources_hook)
)
endif (NOT SuiteSparse_FOUND)
if ((NOT MPI_FOUND) OR (NOT DUNE_ISTL_FOUND))
list (REMOVE_ITEM tests_SOURCES
${PROJECT_SOURCE_DIR}/tests/test_parallel_linearsolver.cpp
)
endif ((NOT MPI_FOUND) OR (NOT DUNE_ISTL_FOUND))
# we are not supposed to include the TinyXML test prog. regardless
list (REMOVE_ITEM opm-core_SOURCES
${PROJECT_SOURCE_DIR}/${opm-core_DIR}/core/utility/parameters/tinyxml/xmltest.cpp
@ -118,3 +121,5 @@ endmacro (install_hook)
# all setup common to the OPM library modules is done here
include (OpmLibMain)

14
CMakeLists_files.cmake
View File

@ -165,9 +165,12 @@ list (APPEND TEST_SOURCE_FILES
tests/test_column_extract.cpp
tests/test_geom2d.cpp
tests/test_linearsolver.cpp
tests/test_parallel_linearsolver.cpp
tests/test_param.cpp
tests/test_blackoilfluid.cpp
tests/test_shadow.cpp
tests/test_equil.cpp
tests/test_regionmapping.cpp
tests/test_units.cpp
tests/test_blackoilstate.cpp
tests/test_parser.cpp
@ -184,6 +187,12 @@ list (APPEND TEST_DATA_FILES
tests/extratestdata.xml
tests/testdata.xml
tests/liveoil.DATA
tests/capillary.DATA
tests/capillary_overlap.DATA
tests/deadfluids.DATA
tests/equil_livegas.DATA
tests/equil_liveoil.DATA
tests/equil_rsvd_and_rvvd.DATA
tests/wetgas.DATA
tests/testBlackoilState1.DATA
tests/testBlackoilState2.DATA
@ -197,6 +206,7 @@ list (APPEND TEST_DATA_FILES
# originally generated with the command:
# find tutorials examples -name '*.c*' -printf '\t%p\n' | sort
list (APPEND EXAMPLE_SOURCE_FILES
examples/compute_initial_state.cpp
examples/compute_tof.cpp
examples/compute_tof_from_files.cpp
examples/import_rewrite.cpp
@ -328,6 +338,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
opm/core/props/satfunc/SaturationPropsInterface.hpp
opm/core/simulator/BlackoilState.hpp
opm/core/simulator/EquilibrationHelpers.hpp
opm/core/simulator/SimulatorCompressibleTwophase.hpp
opm/core/simulator/SimulatorIncompTwophase.hpp
opm/core/simulator/SimulatorOutput.hpp
@ -339,6 +350,8 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/simulator/WellState.hpp
opm/core/simulator/initState.hpp
opm/core/simulator/initState_impl.hpp
opm/core/simulator/initStateEquil.hpp
opm/core/simulator/initStateEquil_impl.hpp
opm/core/tof/DGBasis.hpp
opm/core/tof/TofReorder.hpp
opm/core/tof/TofDiscGalReorder.hpp
@ -373,6 +386,7 @@ list (APPEND PUBLIC_HEADER_FILES
opm/core/utility/NonuniformTableLinear.hpp
opm/core/utility/NullStream.hpp
opm/core/utility/PolynomialUtils.hpp
opm/core/utility/RegionMapping.hpp
opm/core/utility/RootFinders.hpp
opm/core/utility/SparseTable.hpp
opm/core/utility/SparseVector.hpp

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

@ -5,6 +5,8 @@
set (opm-core_CONFIG_VAR
HAVE_ERT
HAVE_SUITESPARSE_UMFPACK_H
HAVE_DUNE_ISTL
HAVE_MPI
)
# dependencies
@ -27,6 +29,8 @@ set (opm-core_DEPS
"TinyXML"
# Ensembles-based Reservoir Tools (ERT)
"ERT"
# Look for MPI support
"MPI"
# DUNE dependency
"dune-common"
"dune-istl"

110
examples/compute_initial_state.cpp Normal file
View File

@ -0,0 +1,110 @@
/*
Copyright 2014 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/>.
*/
#if HAVE_CONFIG_H
#include "config.h"
#endif // HAVE_CONFIG_H
#include <opm/core/grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <boost/filesystem.hpp>
namespace
{
void warnIfUnusedParams(const Opm::parameter::ParameterGroup& param)
{
if (param.anyUnused()) {
std::cout << "-------------------- Unused parameters: --------------------\n";
param.displayUsage();
std::cout << "----------------------------------------------------------------" << std::endl;
}
}
void outputData(const std::string& output_dir,
const std::string& name,
const std::vector<double>& data)
{
std::ostringstream fname;
fname << output_dir << "/" << name;
boost::filesystem::path fpath = fname.str();
try {
create_directories(fpath);
}
catch (...) {
OPM_THROW(std::runtime_error, "Creating directories failed: " << fpath);
}
fname << "/" << "initial.txt";
std::ofstream file(fname.str().c_str());
if (!file) {
OPM_THROW(std::runtime_error, "Failed to open " << fname.str());
}
std::copy(data.begin(), data.end(), std::ostream_iterator<double>(file, "\n"));
}
} // anon namespace
// ----------------- Main program -----------------
int
main(int argc, char** argv)
try
{
using namespace Opm;
// Setup.
parameter::ParameterGroup param(argc, argv, false);
std::cout << "--------------- Reading parameters ---------------" << std::endl;
const std::string deck_filename = param.get<std::string>("deck_filename");
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile(deck_filename);
const double grav = param.getDefault("gravity", unit::gravity);
//EclipseGridParser deck(deck_filename);
GridManager gm(deck);
const UnstructuredGrid& grid = *gm.c_grid();
BlackoilPropertiesFromDeck props(deck, grid, param);
warnIfUnusedParams(param);
// Initialisation.
BlackoilState state;
initStateEquil(grid, props, deck, grav, state);
// Output.
const std::string output_dir = param.getDefault<std::string>("output_dir", "output");
outputData(output_dir, "pressure", state.pressure());
outputData(output_dir, "saturation", state.saturation());
outputData(output_dir, "rs", state.gasoilratio());
outputData(output_dir, "rv", state.rv());
}
catch (const std::exception& e) {
std::cerr << "Program threw an exception: " << e.what() << "\n";
throw;
}

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

@ -25,11 +25,9 @@ struct MultiWriter : public OutputWriter {
MultiWriter (ptr_t writers) : writers_ (std::move (writers)) { }
/// Forward the call to all writers
virtual void writeInit(const SimulatorTimer &timer,
const SimulatorState& reservoirState,
const WellState& wellState) {
virtual void writeInit(const SimulatorTimer &timer) {
for (it_t it = writers_->begin (); it != writers_->end (); ++it) {
(*it)->writeInit (timer, reservoirState, wellState);
(*it)->writeInit (timer);
}
}
@ -48,7 +46,7 @@ private:
/// Psuedo-constructor, can appear in template
template <typename Format> unique_ptr <OutputWriter>
create (const ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const Deck> parser,
std::shared_ptr <const UnstructuredGrid> grid) {
return unique_ptr <OutputWriter> (new Format (params, parser, grid));
}
@ -61,7 +59,7 @@ create (const ParameterGroup& params,
/// to the list below!
typedef map <const char*, unique_ptr <OutputWriter> (*)(
const ParameterGroup&,
std::shared_ptr <EclipseGridParser>,
std::shared_ptr <const Deck>,
std::shared_ptr <const UnstructuredGrid>)> map_t;
map_t FORMATS = {
{ "output_ecl", &create <EclipseWriter> },
@ -71,7 +69,7 @@ map_t FORMATS = {
unique_ptr <OutputWriter>
OutputWriter::create (const ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const Deck> parser,
std::shared_ptr <const UnstructuredGrid> grid) {
// allocate a list which will be filled with writers. this list
// is initially empty (no output).

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

@ -27,7 +27,7 @@ struct UnstructuredGrid;
namespace Opm {
// forward declaration
class EclipseGridParser;
class Deck;
namespace parameter { class ParameterGroup; }
class SimulatorState;
class SimulatorTimer;
@ -43,7 +43,7 @@ class WellState;
* \example
* \code{.cpp}
* ParameterGroup params (argc, argv, false);
* auto parser = std::make_shared <EclipseGridParser> (
* auto parser = std::make_shared <const Deck> (
* params.get <string> ("deck_filename"));
*
* std::unique_ptr <OutputWriter> writer =
@ -67,15 +67,12 @@ public:
virtual ~OutputWriter () { }
/**
* Write the static eclipse data (grid, PVT curves, etc) as well as the
* initial state to disk.
* Write the static data (grid, PVT curves, etc) to disk.
*
* This routine should be called before the first timestep (i.e. when
* timer.currentStepNum () == 0)
*/
virtual void writeInit(const SimulatorTimer &timer,
const SimulatorState& reservoirState,
const WellState& wellState) = 0;
virtual void writeInit(const SimulatorTimer &timer) = 0;
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
@ -88,8 +85,8 @@ public:
* i.e. timer.currentStepNum () > 0.
*/
virtual void writeTimeStep(const SimulatorTimer& timer,
const SimulatorState& reservoirState,
const WellState& wellState) = 0;
const SimulatorState& reservoirState,
const WellState& wellState) = 0;
/*!
* Create a suitable set of output formats based on configuration.
@ -108,7 +105,7 @@ public:
*/
static std::unique_ptr <OutputWriter>
create (const parameter::ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const Deck> parser,
std::shared_ptr <const UnstructuredGrid> grid);
};

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

@ -22,7 +22,6 @@
#include "EclipseWriter.hpp"
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
@ -299,10 +298,6 @@ struct EclipseKeyword : public EclipseHandle <ecl_kw_type> {
static_cast<void> (name);
}
// constructor to keep compatibility with the old eclipse parser
EclipseKeyword (const std::string& name,
const EclipseGridParser& parser);
// GCC 4.4 doesn't generate these constructors for us; provide the
// default implementation explicitly here instead
EclipseKeyword (EclipseKeyword&& rhs)
@ -392,24 +387,6 @@ template <> ecl_type_enum EclipseKeyword<int >::type () { return ECL_INT_TYPE
template <> ecl_type_enum EclipseKeyword<float >::type () { return ECL_FLOAT_TYPE ; }
template <> ecl_type_enum EclipseKeyword<double>::type () { return ECL_DOUBLE_TYPE; }
/// keywords in ERT requires single-precision type, but OPM have them
/// stored as double-precision. this template specialization instantiates
/// a copy function that downcast the data to the required type.
template <>
EclipseKeyword <float>::EclipseKeyword (
const std::string& name,
const EclipseGridParser& parser)
// allocate handle and put in smart pointer base class
: EclipseHandle <ecl_kw_type> (
ecl_kw_alloc (name.c_str(),
// we can safely use the *size* of the original
dataSize (parser.getValue <double> (name), 0, 1),
type ()),
ecl_kw_free) {
const std::vector <double>& data = parser.getValue <double> (name);
copyData (data, &noConversion, 0, 1);
}
/**
* Pointer to memory that holds the name to an Eclipse output file.
*/
@ -435,29 +412,6 @@ private:
}
};
/// Get dimensions of the grid from the parse of the input file
std::vector <int> parserDim (const EclipseGridParser& parser) {
std::vector<int> dim(/* n = */ 3);
// dimensions explicitly given
if (parser.hasField("SPECGRID")) {
dim = parser.getSPECGRID ().dimensions;
}
// dimensions implicitly given by number of deltas
else if (parser.hasField("DXV")) {
assert(parser.hasField("DYV"));
assert(parser.hasField("DZV"));
dim[0] = parser.getFloatingPointValue("DXV").size();
dim[1] = parser.getFloatingPointValue("DYV").size();
dim[2] = parser.getFloatingPointValue("DZV").size();
}
else {
OPM_THROW(std::runtime_error,
"Only decks featureing either the SPECGRID or the D[XYZ]V keywords "
"are currently supported");
}
return dim;
}
/// Get dimensions of the grid from the parse of the input file
std::vector <int> parserDim (Opm::DeckConstPtr newParserDeck) {
std::vector<int> dim(/* n = */ 3);
@ -504,24 +458,6 @@ struct EclipseRestart : public EclipseHandle <ecl_rst_file_type> {
outputStepIdx)),
ecl_rst_file_close) { }
void writeHeader (const SimulatorTimer& timer,
int outputStepIdx,
const PhaseUsage uses,
const EclipseGridParser parser,
const int num_active_cells) {
const std::vector <int> dim = parserDim (parser);
ecl_rst_file_fwrite_header (*this,
outputStepIdx,
timer.currentPosixTime(),
Opm::unit::convert::to (timer.simulationTimeElapsed (),
Opm::unit::day),
dim[0],
dim[1],
dim[2],
num_active_cells,
phaseMask (uses));
}
void writeHeader (const SimulatorTimer& timer,
int outputStepIdx,
const PhaseUsage uses,
@ -568,51 +504,6 @@ private:
* Representation of an Eclipse grid.
*/
struct EclipseGrid : public EclipseHandle <ecl_grid_type> {
/// Create a grid based on the keywords available in input file
static EclipseGrid make (const EclipseGridParser& parser,
int number_of_cells,
const int* cart_dims,
const int* global_cell) {
if (parser.hasField("DXV")) {
// make sure that the DYV and DZV keywords are present if the
// DXV keyword is used in the deck...
assert(parser.hasField("DYV"));
assert(parser.hasField("DZV"));
const auto& dxv = parser.getFloatingPointValue("DXV");
const auto& dyv = parser.getFloatingPointValue("DYV");
const auto& dzv = parser.getFloatingPointValue("DZV");
return EclipseGrid (dxv, dyv, dzv);
}
else if (parser.hasField("ZCORN")) {
struct grdecl g = parser.get_grdecl ();
auto coordData = parser.getFloatingPointValue(COORD_KW);
auto zcornData = parser.getFloatingPointValue(ZCORN_KW);
EclipseKeyword<float> coord_kw (COORD_KW, coordData);
EclipseKeyword<float> zcorn_kw (ZCORN_KW, zcornData);
// get the actually active cells, after processing
std::vector <int> actnum;
getActiveCells_(number_of_cells, cart_dims, global_cell, actnum);
EclipseKeyword<int> actnum_kw (ACTNUM_KW, actnum);
EclipseKeyword<float> mapaxes_kw (MAPAXES_KW);
if (g.mapaxes) {
auto mapaxesData = parser.getFloatingPointValue(MAPAXES_KW);
mapaxes_kw = std::move (EclipseKeyword<float> (MAPAXES_KW, mapaxesData));
}
return EclipseGrid (cart_dims, zcorn_kw, coord_kw, actnum_kw, mapaxes_kw);
}
else {
OPM_THROW(std::runtime_error,
"Can't create an ERT grid (no supported keywords found in deck)");
}
}
/// Create a grid based on the keywords available in input file
static EclipseGrid make (Opm::DeckConstPtr newParserDeck,
int number_of_cells,
@ -743,18 +634,6 @@ struct EclipseInit : public EclipseHandle <fortio_type> {
return EclipseInit (initFileName, fmt_file);
}
void writeHeader (const EclipseGrid& grid,
const SimulatorTimer& timer,
const EclipseGridParser& parser,
const PhaseUsage uses) {
EclipseKeyword<float> poro (PORO_KW, parser);
ecl_init_file_fwrite_header (*this,
grid,
poro,
phaseMask (uses),
timer.currentPosixTime());
}
void writeHeader (int number_of_cells,
const int* cart_dims,
const int* global_cell,
@ -776,15 +655,6 @@ struct EclipseInit : public EclipseHandle <fortio_type> {
timer.currentPosixTime ());
}
void writeKeyword (const std::string& keywordName,
const EclipseGridParser& parser,
double (* const transf)(const double&)) {
auto data = parser.getValue <double> (keywordName);
convertUnit_(data, transf);
EclipseKeyword <float> kw (keywordName, data);
ecl_kw_fwrite (kw, *this);
}
void writeKeyword (const std::string& keywordName, const std::vector<double> &data)
{
EclipseKeyword <float> kw (keywordName, data);
@ -824,14 +694,6 @@ template struct EclipseHandle<ecl_sum_tstep_struct>;
struct EclipseWellReport;
struct EclipseSummary : public EclipseHandle <ecl_sum_type> {
EclipseSummary (const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
const EclipseGridParser parser)
: EclipseHandle <ecl_sum_type> (
alloc_writer (outputDir, baseName, timer, parser),
ecl_sum_free) { }
EclipseSummary (const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
@ -855,10 +717,6 @@ struct EclipseSummary : public EclipseHandle <ecl_sum_type> {
ecl_sum_fwrite (*this);
}
// add rate variables for each of the well in the input file
void addWells (const EclipseGridParser& parser,
const PhaseUsage& uses);
// add rate variables for each of the well in the input file
void addWells (Opm::DeckConstPtr newParserDeck,
const PhaseUsage& uses);
@ -889,27 +747,6 @@ private:
return std::unique_ptr <EclipseTimeStep> (tstep);
}
/// Helper routine that lets us use local variables to hold
/// intermediate results while filling out the allocations function's
/// argument list.
static ecl_sum_type* alloc_writer (const std::string& outputDir,
const std::string& baseName,
const SimulatorTimer& timer,
const EclipseGridParser& parser) {
boost::filesystem::path casePath (outputDir);
casePath /= boost::to_upper_copy (baseName);
const std::vector <int> dim = parserDim (parser);
return ecl_sum_alloc_writer (casePath.string ().c_str (),
false, /* formatted */
true, /* unified */
":", /* join string */
timer.simulationTimeElapsed (),
dim[0],
dim[1],
dim[2]);
}
/// Helper routine that lets us use local variables to hold
/// intermediate results while filling out the allocations function's
/// argument list.
@ -938,29 +775,6 @@ private:
*/
struct EclipseWellReport : public EclipseHandle <smspec_node_type> {
protected:
EclipseWellReport (const EclipseSummary& summary, /* section to add to */
const EclipseGridParser& parser, /* well names */
int whichWell, /* index of well line */
PhaseUsage uses, /* phases present */
BlackoilPhases::PhaseIndex phase, /* oil, water or gas */
WellType type, /* prod. or inj. */
char aggregation, /* rate or total or BHP */
std::string unit)
: EclipseHandle <smspec_node_type> (
ecl_sum_add_var (summary,
varName (phase,
type,
aggregation).c_str (),
wellName (parser, whichWell).c_str (),
/* num = */ 0,
unit.c_str(),
/* defaultValue = */ 0.))
// save these for when we update the value in a timestep
, index_ (whichWell * uses.num_phases + uses.phase_pos [phase])
// producers can be seen as negative injectors
, sign_ (type == INJECTOR ? +1. : -1.) { }
EclipseWellReport (const EclipseSummary& summary, /* section to add to */
Opm::DeckConstPtr newParserDeck, /* well names */
int whichWell, /* index of well line */
@ -1002,12 +816,6 @@ private:
/// natural sign of the rate
const double sign_;
/// Get the name associated with this well
std::string wellName (const EclipseGridParser& parser,
int whichWell) {
return parser.getWELSPECS().welspecs[whichWell].name_;
}
/// Get the name associated with this well
std::string wellName (Opm::DeckConstPtr newParserDeck,
int whichWell)
@ -1063,21 +871,6 @@ protected:
/// Monitors the rate given by a well.
struct EclipseWellRate : public EclipseWellReport {
EclipseWellRate (const EclipseSummary& summary,
const EclipseGridParser& parser,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
parser,
whichWell,
uses,
phase,
type,
'R',
"SM3/DAY" /* surf. cub. m. per day */ ) { }
EclipseWellRate (const EclipseSummary& summary,
Opm::DeckConstPtr newParserDeck,
int whichWell,
@ -1102,24 +895,6 @@ struct EclipseWellRate : public EclipseWellReport {
/// Monitors the total production in a well.
struct EclipseWellTotal : public EclipseWellReport {
EclipseWellTotal (const EclipseSummary& summary,
const EclipseGridParser& parser,
int whichWell,
PhaseUsage uses,
BlackoilPhases::PhaseIndex phase,
WellType type)
: EclipseWellReport (summary,
parser,
whichWell,
uses,
phase,
type,
'T',
"SM3" /* surface cubic meter */ )
// nothing produced when the reporting starts
, total_ (0.) { }
EclipseWellTotal (const EclipseSummary& summary,
Opm::DeckConstPtr newParserDeck,
int whichWell,
@ -1162,22 +937,6 @@ private:
/// 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")
{ }
EclipseWellBhp (const EclipseSummary& summary,
Opm::DeckConstPtr newParserDeck,
int whichWell,
@ -1212,75 +971,15 @@ EclipseSummary::writeTimeStep (const SimulatorTimer& timer,
const double value = (*v)->update (timer, wellState);
ecl_sum_tstep_iset(*tstep, *(*v).get (), value);
}
// write the summary file to disk
ecl_sum_fwrite(*this);
}
/// Supported well types. Enumeration doesn't let us get all the members,
/// so we must have an explicit array.
static WellType WELL_TYPES[] = { INJECTOR, PRODUCER };
inline void
EclipseSummary::addWells (const EclipseGridParser& parser,
const PhaseUsage& uses)
{
// TODO: Only create report variables that are requested with keywords
// (e.g. "WOPR") in the input files, and only for those wells that are
// mentioned in those keywords
const int numWells = parser.getWELSPECS().welspecs.size();
for (int phaseCounter = 0;
phaseCounter != BlackoilPhases::MaxNumPhases;
++phaseCounter) {
const BlackoilPhases::PhaseIndex phase =
static_cast <BlackoilPhases::PhaseIndex> (phaseCounter);
// don't bother with reporting for phases that aren't there
if (!uses.phase_used [phaseCounter]) {
continue;
}
for (size_t typeIndex = 0;
typeIndex < sizeof (WELL_TYPES) / sizeof (WELL_TYPES[0]);
++typeIndex) {
const WellType type = WELL_TYPES[typeIndex];
for (int whichWell = 0; whichWell != numWells; ++whichWell) {
// W{O,G,W}{I,P}R
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellRate (*this,
parser,
whichWell,
uses,
phase,
type)));
// W{O,G,W}{I,P}T
add (std::unique_ptr <EclipseWellReport> (
new EclipseWellTotal (*this,
parser,
whichWell,
uses,
phase,
type)));
}
}
}
// 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])));
}
}
inline void
EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck,
const PhaseUsage& uses) {
@ -1347,171 +1046,52 @@ EclipseSummary::addWells (Opm::DeckConstPtr newParserDeck,
namespace Opm {
void EclipseWriter::writeInit(const SimulatorTimer &timer,
const SimulatorState& reservoirState,
const WellState& wellState)
void EclipseWriter::writeInit(const SimulatorTimer &timer)
{
// if we don't want to write anything, this method becomes a
// no-op...
if (!enableOutput_) {
return;
}
if (newParserDeck_) {
/* Grid files */
EclipseGrid eclGrid = EclipseGrid::make (newParserDeck_, number_of_cells_,
cart_dims_, global_cell_);
eclGrid.write (outputDir_, baseName_, /*stepIdx=*/0);
/* Grid files */
EclipseGrid eclGrid = EclipseGrid::make (newParserDeck_, number_of_cells_,
cart_dims_, global_cell_);
eclGrid.write (outputDir_, baseName_, /*stepIdx=*/0);
EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0);
fortio.writeHeader (number_of_cells_,
cart_dims_,
global_cell_,
timer,
newParserDeck_,
uses_);
EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0);
fortio.writeHeader (number_of_cells_,
cart_dims_,
global_cell_,
timer,
newParserDeck_,
uses_);
if (newParserDeck_->hasKeyword("PERM")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERM"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERM", data);
}
if (newParserDeck_->hasKeyword("PERMX")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMX"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMX", data);
}
if (newParserDeck_->hasKeyword("PERMY")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMY"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMY", data);
}
if (newParserDeck_->hasKeyword("PERMZ")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMZ"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMZ", data);
}
/* Initial solution (pressure and saturation) */
writeSolution_(timer, reservoirState);
/* 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, newParserDeck_));
summary_->addWells (newParserDeck_, uses_);
}
else {
/* Grid files */
EclipseGrid ecl_grid = EclipseGrid::make (*parser_, number_of_cells_,
cart_dims_, global_cell_);
ecl_grid.write (outputDir_, baseName_, /*stepIdx=*/0);
EclipseInit fortio = EclipseInit::make (outputDir_, baseName_, /*stepIdx=*/0);
fortio.writeHeader (ecl_grid,
timer,
*parser_,
uses_);
fortio.writeKeyword ("PERMX", *parser_, &toMilliDarcy);
fortio.writeKeyword ("PERMY", *parser_, &toMilliDarcy);
fortio.writeKeyword ("PERMZ", *parser_, &toMilliDarcy);
/* Initial solution (pressure and saturation) */
writeSolution_(timer, reservoirState);
/* 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::activeToGlobalCellData_(std::vector<double> &globalCellsBuf,
const std::vector<double> &activeCellsBuf,
const std::vector<double> &inactiveCellsBuf) const
{
globalCellsBuf = inactiveCellsBuf;
// overwrite the values of active cells
for (int activeCellIdx = 0;
activeCellIdx < number_of_cells_;
++activeCellIdx)
{
int globalCellIdx = global_cell_[activeCellIdx];
globalCellsBuf[globalCellIdx] = activeCellsBuf[activeCellIdx];
}
}
void EclipseWriter::writeSolution_(const SimulatorTimer& timer,
const SimulatorState& reservoirState)
{
if (newParserDeck_) {
// start writing to files
EclipseRestart rst(outputDir_, baseName_, timer, outputTimeStepIdx_);
rst.writeHeader (timer, outputTimeStepIdx_, uses_, newParserDeck_, reservoirState.pressure().size ());
EclipseSolution sol (rst);
// write out the pressure of the reference phase (whatever
// phase that is...). this is not the most performant solution
// thinkable, but this is also not in the most performance
// critical code path!
std::vector<double> tmp = reservoirState.pressure();
convertUnit_(tmp, toBar);
sol.add(EclipseKeyword<float>("PRESSURE", tmp));
for (int phase = 0; phase != BlackoilPhases::MaxNumPhases; ++phase) {
// Eclipse never writes the oil saturation, so all post-processors
// must calculate this from the other saturations anyway
if (phase == BlackoilPhases::PhaseIndex::Liquid) {
continue;
}
if (uses_.phase_used [phase]) {
tmp = reservoirState.saturation();
extractFromStripedData_(tmp,
/*offset=*/uses_.phase_pos[phase],
/*stride=*/uses_.num_phases);
sol.add(EclipseKeyword<float>(SAT_NAMES[phase], tmp));
}
}
}
else {
// start writing to files
EclipseRestart rst (outputDir_,
baseName_,
timer,
outputTimeStepIdx_);
rst.writeHeader (timer,
outputTimeStepIdx_,
uses_,
*parser_,
reservoirState.pressure ().size ());
EclipseSolution sol (rst);
// write pressure and saturation fields (same as DataMap holds)
// convert the pressures from Pascals to bar because Eclipse
// seems to write bars
auto data = reservoirState.pressure();
convertUnit_(data, toBar);
sol.add(EclipseKeyword<float>("PRESSURE", data));
for (int phase = 0; phase != BlackoilPhases::MaxNumPhases; ++phase) {
// Eclipse never writes the oil saturation, so all post-processors
// must calculate this from the other saturations anyway
if (phase == BlackoilPhases::PhaseIndex::Liquid) {
continue;
}
if (uses_.phase_used [phase]) {
auto tmp = reservoirState.saturation();
extractFromStripedData_(tmp,
/*offset=*/uses_.phase_pos[phase],
/*stride=*/uses_.num_phases);
sol.add (EclipseKeyword<float>(SAT_NAMES[phase], tmp));
}
}
if (newParserDeck_->hasKeyword("PERM")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERM"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERM", data);
}
++outputTimeStepIdx_;
if (newParserDeck_->hasKeyword("PERMX")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMX"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMX", data);
}
if (newParserDeck_->hasKeyword("PERMY")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMY"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMY", data);
}
if (newParserDeck_->hasKeyword("PERMZ")) {
auto data = getAllSiDoubles_(newParserDeck_->getKeyword("PERMZ"));
convertUnit_(data, toMilliDarcy);
fortio.writeKeyword ("PERMZ", data);
}
/* 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, newParserDeck_));
summary_->addWells (newParserDeck_, uses_);
}
void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
@ -1525,12 +1105,38 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
}
// respected the output_interval parameter
if (timer.currentStepNum() % outputInterval_ != 0) {
if (outputTimeStepIdx_ % outputInterval_ != 0) {
return;
}
/* Field variables (pressure, saturation) */
writeSolution_(timer, reservoirState);
// start writing to files
EclipseRestart rst(outputDir_, baseName_, timer, outputTimeStepIdx_);
rst.writeHeader (timer, outputTimeStepIdx_, uses_, newParserDeck_, reservoirState.pressure().size ());
EclipseSolution sol (rst);
// write out the pressure of the reference phase (whatever
// phase that is...). this is not the most performant solution
// thinkable, but this is also not in the most performance
// critical code path!
std::vector<double> tmp = reservoirState.pressure();
convertUnit_(tmp, toBar);
sol.add(EclipseKeyword<float>("PRESSURE", tmp));
for (int phase = 0; phase != BlackoilPhases::MaxNumPhases; ++phase) {
// Eclipse never writes the oil saturation, so all post-processors
// must calculate this from the other saturations anyway
if (phase == BlackoilPhases::PhaseIndex::Liquid) {
continue;
}
if (uses_.phase_used [phase]) {
tmp = reservoirState.saturation();
extractFromStripedData_(tmp,
/*offset=*/uses_.phase_pos[phase],
/*stride=*/uses_.num_phases);
sol.add(EclipseKeyword<float>(SAT_NAMES[phase], tmp));
}
}
/* Summary variables (well reporting) */
// TODO: instead of writing the header (smspec) every time, it should
@ -1547,14 +1153,14 @@ void EclipseWriter::writeTimeStep(const SimulatorTimer& timer,
// will contain data from the whole simulation, instead of just
// the last step.
summary_->writeTimeStep(timer, wellState);
++outputTimeStepIdx_;
}
#else
namespace Opm {
void EclipseWriter::writeInit(const SimulatorTimer&,
const SimulatorState&,
const WellState&)
void EclipseWriter::writeInit(const SimulatorTimer&)
{
// if we don't want to write anything, this method becomes a
// no-op...
@ -1585,93 +1191,31 @@ void EclipseWriter::writeTimeStep(
EclipseWriter::EclipseWriter (
const ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid)
: parser_ (parser)
, newParserDeck_()
, number_of_cells_(grid->number_of_cells)
, dimensions_(grid->dimensions)
, cart_dims_(grid->cartdims)
, global_cell_(grid->global_cell)
, uses_ (phaseUsageFromDeck (*parser))
{
init(params);
}
EclipseWriter::EclipseWriter (
const ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell, const int* cart_dims,
int dimensions)
: parser_ (parser)
: newParserDeck_ (newParserDeck)
, number_of_cells_(number_of_cells)
, dimensions_(dimensions)
, cart_dims_(cart_dims)
, global_cell_(global_cell)
, uses_ (phaseUsageFromDeck (*parser)) {
, uses_ (phaseUsageFromDeck (newParserDeck_)) {
init(params);
}
void EclipseWriter::init(const ParameterGroup& params)
{
// set the index of the first time step written to 0...
outputTimeStepIdx_ = 0;
// get the base name from the name of the deck
using boost::filesystem::path;
path deck (params.get <std::string> ("deck_filename"));
if (boost::to_upper_copy (path (deck.extension ()).string ()) == ".DATA") {
baseName_ = path (deck.stem ()).string ();
}
else {
baseName_ = path (deck.filename ()).string ();
}
// make uppercase of everything (or otherwise we'll get uppercase
// of some of the files (.SMSPEC, .UNSMRY) and not others
baseName_ = boost::to_upper_copy (baseName_);
// retrieve the value of the "output" parameter
enableOutput_ = params.getDefault<bool>("output", /*defaultValue=*/true);
// retrieve the interval at which something should get written to
// disk (once every N timesteps)
outputInterval_ = params.getDefault<int>("output_interval", /*defaultValue=*/1);
// store in current directory if not explicitly set
outputDir_ = params.getDefault<std::string>("output_dir", ".");
// set the index of the first time step written to 0...
outputTimeStepIdx_ = 0;
if (enableOutput_) {
// make sure that the output directory exists, if not try to create it
if (!boost::filesystem::exists(outputDir_)) {
std::cout << "Trying to create directory \"" << outputDir_ << "\" for the simulation output\n";
boost::filesystem::create_directories(outputDir_);
}
if (!boost::filesystem::is_directory(outputDir_)) {
OPM_THROW(std::runtime_error,
"The path specified as output directory '" << outputDir_
<< "' is not a directory");
}
}
}
EclipseWriter::EclipseWriter (
const ParameterGroup& params,
Opm::DeckConstPtr newParserDeck,
std::shared_ptr <const UnstructuredGrid> grid)
: parser_ ()
, newParserDeck_(newParserDeck)
std::shared_ptr<const UnstructuredGrid> grid)
: newParserDeck_ (newParserDeck)
, number_of_cells_(grid->number_of_cells)
, dimensions_(grid->dimensions)
, cart_dims_(grid->cartdims)
, global_cell_(grid->global_cell)
, uses_ (phaseUsageFromDeck (newParserDeck))
, uses_ (phaseUsageFromDeck (newParserDeck_)) {
init(params);
}
void EclipseWriter::init(const ParameterGroup& params)
{
// get the base name from the name of the deck
using boost::filesystem::path;

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

@ -36,7 +36,6 @@ struct EclipseSummary;
namespace Opm {
// forward declarations
class EclipseGridParser;
class SimulatorState;
class SimulatorTimer;
class WellState;
@ -55,19 +54,6 @@ namespace parameter { class ParameterGroup; }
class EclipseWriter : public OutputWriter
{
public:
/*!
* \brief Sets the common attributes required to write eclipse
* binary files using ERT.
*/
EclipseWriter(const parameter::ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const UnstructuredGrid> grid);
EclipseWriter (const parameter::ParameterGroup& params,
std::shared_ptr <const EclipseGridParser> parser,
int number_of_cells, const int* global_cell, const int* cart_dims,
int dimension);
/*!
* \brief Sets the common attributes required to write eclipse
* binary files using ERT.
@ -94,16 +80,20 @@ public:
virtual ~EclipseWriter ();
/**
* Write the static eclipse data (grid, PVT curves, etc) as well as the
* initial state to disk.
* Write the static eclipse data (grid, PVT curves, etc) to disk.
*/
virtual void writeInit(const SimulatorTimer &timer,
const SimulatorState& reservoirState,
const WellState& wellState);
virtual void writeInit(const SimulatorTimer &timer);
/*!
* \brief Write a blackoil reservoir state to disk for later inspection with
* visualization tools like ResInsight
* \brief Write a reservoir state and summary information to disk.
*
*
* The reservoir state can be inspected with visualization tools like
* ResInsight.
*
* The summary information can then be visualized using tools from
* ERT or ECLIPSE. Note that calling this method is only
* meaningful after the first time step has been completed.
*
* \param[in] reservoirState The thermodynamic state of the reservoir
* \param[in] wellState The production/injection data for all wells
@ -113,7 +103,6 @@ public:
const WellState& wellState);
private:
std::shared_ptr <const EclipseGridParser> parser_;
Opm::DeckConstPtr newParserDeck_;
int number_of_cells_;
int dimensions_;
@ -127,14 +116,6 @@ private:
PhaseUsage uses_; // active phases in the input deck
std::shared_ptr <EclipseSummary> summary_;
void activeToGlobalCellData_(std::vector<double> &globalCellsBuf,
const std::vector<double> &activeCellsBuf,
const std::vector<double> &inactiveCellsBuf) const;
/// Write solution field variables (pressure and saturation)
void writeSolution_(const SimulatorTimer& timer,
const SimulatorState& reservoirState);
void init(const parameter::ParameterGroup& params);
};
} // namespace Opm

5
opm/core/linalg/LinearSolverFactory.cpp
View File

@ -97,9 +97,10 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const
double* solution,
const boost::any& add) const
{
return solver_->solve(size, nonzeros, ia, ja, sa, rhs, solution);
return solver_->solve(size, nonzeros, ia, ja, sa, rhs, solution, add);
}
void LinearSolverFactory::setTolerance(const double tol)

3
opm/core/linalg/LinearSolverFactory.hpp
View File

@ -74,7 +74,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const;
double* solution,
const boost::any& add=boost::any()) const;
/// Set tolerance for the linear solver.
/// \param[in] tol tolerance value

4
opm/core/linalg/LinearSolverInterface.hpp
View File

@ -20,6 +20,7 @@
#ifndef OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED
#define OPM_LINEARSOLVERINTERFACE_HEADER_INCLUDED
#include<boost/any.hpp>
struct CSRMatrix;
@ -69,7 +70,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const = 0;
double* solution,
const boost::any& add=boost::any()) const = 0;
/// Set tolerance for the linear solver.
/// \param[in] tol tolerance value

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

@ -23,6 +23,8 @@
#endif
#include <opm/core/linalg/LinearSolverIstl.hpp>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#include <opm/core/utility/ErrorMacros.hpp>
// Silence compatibility warning from DUNE headers since we don't use
// the deprecated member anyway (in this compilation unit)
@ -35,7 +37,9 @@
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/operators.hh>
#include <dune/istl/io.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/istl/preconditioners.hh>
#include <dune/istl/schwarz.hh>
#include <dune/istl/solvers.hh>
#include <dune/istl/paamg/amg.hh>
#include <dune/istl/paamg/kamg.hh>
@ -47,7 +51,7 @@
#include <stdexcept>
#include <iostream>
#include <type_traits>
namespace Opm
{
@ -134,7 +138,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const
double* solution,
const boost::any& comm) const
{
// Build Istl structures from input.
// System matrix
@ -155,10 +160,26 @@ namespace Opm
if (maxit == 0) {
maxit = 5000;
}
Dune::SeqScalarProduct<Vector> sp;
Dune::Amg::SequentialInformation comm;
Operator opA(A);
return solveSystem(opA, solution, rhs, sp, comm, maxit);
#if HAVE_MPI
if(comm.type()==typeid(ParallelISTLInformation))
{
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
const ParallelISTLInformation& info = boost::any_cast<const ParallelISTLInformation&>(comm);
Comm istlComm(info.communicator());
info.copyValuesTo(istlComm.indexSet(), istlComm.remoteIndices());
Dune::OverlappingSchwarzOperator<Mat,Vector,Vector, Comm>
opA(A, istlComm);
Dune::OverlappingSchwarzScalarProduct<Vector,Comm> sp(istlComm);
return solveSystem(opA, solution, rhs, sp, istlComm, maxit);
}
else
#endif
{
Dune::SeqScalarProduct<Vector> sp;
Dune::Amg::SequentialInformation comm;
Operator opA(A);
return solveSystem(opA, solution, rhs, sp, comm, maxit);
}
}
template<class O, class S, class C>
@ -204,6 +225,14 @@ namespace Opm
break;
case FastAMG:
#if defined(HAS_DUNE_FAST_AMG) || DUNE_VERSION_NEWER(DUNE_ISTL, 2, 3)
#if HAVE_MPI
if(std::is_same<C,Dune::OwnerOverlapCopyCommunication<int,int> >::value)
{
OPM_THROW(std::runtime_error, "Trying to use sequential FastAMG solver for a parallel problem!");
}
#endif // HAVE_MPI
res = solveFastAMG(opA, x, b, sp, comm, linsolver_residual_tolerance_, maxit, linsolver_verbosity_,
linsolver_prolongate_factor_);
#else
@ -234,28 +263,47 @@ namespace Opm
return linsolver_residual_tolerance_;
}
namespace
{
template<class P, class O, class C>
struct SmootherChooser
{
typedef P Type;
};
#if HAVE_MPI
template<class P, class O>
struct SmootherChooser<P, O, Dune::OwnerOverlapCopyCommunication<int,int> >
{
typedef Dune::OwnerOverlapCopyCommunication<int,int> Comm;
typedef Dune::BlockPreconditioner<typename O::domain_type, typename O::range_type,
Comm, P>
Type;
};
#endif
template<class P, class O, class C>
struct PreconditionerTraits
{
typedef std::shared_ptr<P> PointerType;
typedef typename SmootherChooser<P,O,C>::Type SmootherType;
typedef std::shared_ptr<SmootherType> PointerType;
};
template<class P, class O, class C>
typename PreconditionerTraits<P,O>::PointerType
typename PreconditionerTraits<P,O,C>::PointerType
makePreconditioner(O& opA, double relax, const C& comm, int iterations=1)
{
typename Dune::Amg::SmootherTraits<P>::Arguments args;
typename Dune::Amg::ConstructionTraits<P>::Arguments cargs;
typedef typename SmootherChooser<P,O,C>::Type SmootherType;
typedef typename PreconditionerTraits<P,O,C>::PointerType PointerType;
typename Dune::Amg::SmootherTraits<SmootherType>::Arguments args;
typename Dune::Amg::ConstructionTraits<SmootherType>::Arguments cargs;
cargs.setMatrix(opA.getmat());
args.iterations=iterations;
args.relaxationFactor=relax;
cargs.setArgs(args);
cargs.setComm(comm);
return std::shared_ptr<P>(Dune::Amg::ConstructionTraits<P>::construct(cargs));
return PointerType(Dune::Amg::ConstructionTraits<SmootherType>::construct(cargs));
}
template<class O, class S, class C>
@ -323,19 +371,20 @@ namespace Opm
#endif
#if SMOOTHER_ILU
typedef Dune::SeqILU0<Mat,Vector,Vector> Smoother;
typedef Dune::SeqILU0<Mat,Vector,Vector> SeqSmoother;
#else
typedef Dune::SeqSOR<Mat,Vector,Vector> Smoother;
typedef Dune::SeqSOR<Mat,Vector,Vector> SeqSmoother;
#endif
typedef typename SmootherChooser<SeqSmoother, O, C>::Type Smoother;
typedef Dune::Amg::CoarsenCriterion<CriterionBase> Criterion;
typedef Dune::Amg::AMG<Operator,Vector,Smoother> Precond;
typedef Dune::Amg::AMG<O,Vector,Smoother,C> Precond;
// Construct preconditioner.
Criterion criterion;
Precond::SmootherArgs smootherArgs;
typename Precond::SmootherArgs smootherArgs;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
linsolver_smooth_steps);
Precond precond(opA, criterion, smootherArgs);
Precond precond(opA, criterion, smootherArgs, comm);
// Construct linear solver.
Dune::CGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
@ -360,6 +409,7 @@ namespace Opm
double linsolver_prolongate_factor, int linsolver_smooth_steps)
{
// Solve with AMG solver.
Dune::MatrixAdapter<typename O::matrix_type,Vector,Vector> sOpA(opA.getmat());
#if FIRST_DIAGONAL
typedef Dune::Amg::FirstDiagonal CouplingMetric;
@ -386,10 +436,10 @@ namespace Opm
Criterion criterion;
setUpCriterion(criterion, linsolver_prolongate_factor, verbosity,
linsolver_smooth_steps);
Precond precond(opA, criterion, smootherArgs);
Precond precond(sOpA, criterion, smootherArgs);
// Construct linear solver.
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
Dune::GeneralizedPCGSolver<Vector> linsolve(sOpA, precond, tolerance, maxit, verbosity);
// Solve system.
Dune::InverseOperatorResult result;
@ -409,6 +459,8 @@ namespace Opm
double linsolver_prolongate_factor)
{
// Solve with AMG solver.
typedef Dune::MatrixAdapter<typename O::matrix_type, Vector, Vector> Operator;
Operator sOpA(opA.getmat());
#if FIRST_DIAGONAL
typedef Dune::Amg::FirstDiagonal CouplingMetric;
@ -434,10 +486,10 @@ namespace Opm
parms.setNoPreSmoothSteps(smooth_steps);
parms.setNoPostSmoothSteps(smooth_steps);
parms.setProlongationDampingFactor(linsolver_prolongate_factor);
Precond precond(opA, criterion, parms);
Precond precond(sOpA, criterion, parms);
// Construct linear solver.
Dune::GeneralizedPCGSolver<Vector> linsolve(opA, sp, precond, tolerance, maxit, verbosity);
Dune::GeneralizedPCGSolver<Vector> linsolve(sOpA, precond, tolerance, maxit, verbosity);
// Solve system.
Dune::InverseOperatorResult result;

6
opm/core/linalg/LinearSolverIstl.hpp
View File

@ -24,12 +24,11 @@
#include <opm/core/linalg/LinearSolverInterface.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <string>
#include <boost/any.hpp>
namespace Opm
{
/// Concrete class encapsulating some dune-istl linear solvers.
class LinearSolverIstl : public LinearSolverInterface
{
@ -75,7 +74,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const;
double* solution,
const boost::any& comm=boost::any()) const;
/// Set tolerance for the residual in dune istl linear solver.
/// \param[in] tol tolerance value

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

@ -46,7 +46,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const
double* solution,
const boost::any&) const
{
CSRMatrix A = {
(size_t)size,

3
opm/core/linalg/LinearSolverUmfpack.hpp
View File

@ -55,7 +55,8 @@ namespace Opm
const int* ja,
const double* sa,
const double* rhs,
double* solution) const;
double* solution,
const boost::any& add=boost::any()) const;
/// Set tolerance for the linear solver.
/// \param[in] tol tolerance value

160
opm/core/linalg/ParallelIstlInformation.hpp Normal file
View File

@ -0,0 +1,160 @@
/*
Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
#define OPM_PARALLELISTLINFORMTION_HEADER_INCLUDED
#if HAVE_MPI && HAVE_DUNE_ISTL
#include "mpi.h"
#include <dune/istl/owneroverlapcopy.hh>
#include <dune/common/parallel/interface.hh>
#include <dune/common/parallel/communicator.hh>
#include <dune/common/enumset.hh>
#include<algorithm>
namespace Opm
{
/// \brief Class that encapsulates the parallelization information needed by the
/// ISTL solvers.
class ParallelISTLInformation
{
public:
/// \brief The type of the parallel index set used.
typedef Dune::OwnerOverlapCopyCommunication<int, int>::ParallelIndexSet ParallelIndexSet;
/// \brief The type of the remote indices information used.
typedef Dune::OwnerOverlapCopyCommunication<int, int>::RemoteIndices RemoteIndices;
/// \brief Constructs an empty parallel information object using MPI_COMM_WORLD
ParallelISTLInformation()
: indexSet_(new ParallelIndexSet),
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, MPI_COMM_WORLD)),
communicator_(MPI_COMM_WORLD)
{}
/// \brief Constructs an empty parallel information object using a communicator.
/// \param communicator The communicator to use.
ParallelISTLInformation(MPI_Comm communicator)
: indexSet_(new ParallelIndexSet),
remoteIndices_(new RemoteIndices(*indexSet_, *indexSet_, communicator)),
communicator_(communicator)
{}
/// \brief Constructs a parallel information object from the specified information.
/// \param indexSet The parallel index set to use.
/// \param remoteIndices The remote indices information to use.
/// \param communicator The communicator to use.
ParallelISTLInformation(const std::shared_ptr<ParallelIndexSet>& indexSet,
const std::shared_ptr<RemoteIndices>& remoteIndices,
MPI_Comm communicator)
: indexSet_(indexSet), remoteIndices_(remoteIndices), communicator_(communicator)
{}
/// \brief Copy constructor.
///
/// The information will be shared by the the two objects.
ParallelISTLInformation(const ParallelISTLInformation& other)
: indexSet_(other.indexSet_), remoteIndices_(other.remoteIndices_),
communicator_(other.communicator_)
{}
/// \brief Get a pointer to the underlying index set.
std::shared_ptr<ParallelIndexSet> indexSet() const
{
return indexSet_;
}
/// \brief Get a pointer to the remote indices information.
std::shared_ptr<RemoteIndices> remoteIndices() const
{
return remoteIndices_;
}
/// \brief Get the MPI communicator that we use.
MPI_Comm communicator() const
{
return communicator_;
}
/// \brief Copy the information stored to the specified objects.
/// \param[out] indexSet The object to store the index set in.
/// \param[out] remoteIndices The object to store the remote indices information in.
void copyValuesTo(ParallelIndexSet& indexSet, RemoteIndices& remoteIndices) const
{
indexSet.beginResize();
IndexSetInserter<ParallelIndexSet> inserter(indexSet);
std::for_each(indexSet_->begin(), indexSet_->end(), inserter);
indexSet.endResize();
remoteIndices.rebuild<false>();
}
/// \brief Communcate the dofs owned by us to the other process.
///
/// Afterwards all associated dofs will contain the same data.
template<class T>
void copyOwnerToAll (const T& source, T& dest) const
{
typedef Dune::Combine<Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::owner>,Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::overlap>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> OwnerOverlapSet;
typedef Dune::Combine<OwnerOverlapSet, Dune::EnumItem<Dune::OwnerOverlapCopyAttributeSet::AttributeSet,Dune::OwnerOverlapCopyAttributeSet::copy>,Dune::OwnerOverlapCopyAttributeSet::AttributeSet> AllSet;
OwnerOverlapSet sourceFlags;
AllSet destFlags;
Dune::Interface interface(communicator_);
if(!remoteIndices_->isSynced())
remoteIndices_->rebuild<false>();
interface.build(*remoteIndices_,sourceFlags,destFlags);
Dune::BufferedCommunicator communicator;
communicator.template build<T>(interface);
communicator.template forward<CopyGatherScatter<T> >(source,dest);
communicator.free();
}
private:
/** \brief gather/scatter callback for communcation */
template<typename T>
struct CopyGatherScatter
{
typedef typename Dune::CommPolicy<T>::IndexedType V;
static V gather(const T& a, std::size_t i)
{
return a[i];
}
static void scatter(T& a, V v, std::size_t i)
{
a[i] = v;
}
};
template<class T>
class IndexSetInserter
{
public:
typedef T ParallelIndexSet;
IndexSetInserter(ParallelIndexSet& indexSet)
: indexSet_(&indexSet)
{}
void operator()(const typename ParallelIndexSet::IndexPair& pair)
{
indexSet_->add(pair.global(), pair.local());
}
private:
ParallelIndexSet* indexSet_;
};
std::shared_ptr<ParallelIndexSet> indexSet_;
std::shared_ptr<RemoteIndices> remoteIndices_;
MPI_Comm communicator_;
};
} // end namespace Opm
#endif
#endif

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

@ -34,6 +34,7 @@
#include <opm/core/utility/linearInterpolation.hpp>
#include <algorithm>
#include <iostream>
namespace Opm
{
@ -99,7 +100,15 @@ namespace Opm
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();
undersat_gas_tables_[i][2] = undersatTable.getGasViscosityColumn();
}
// Bg -> 1/Bg
for (int i=0; i<sz; ++i) {
saturated_gas_table_[1][i] = 1.0/saturated_gas_table_[1][i];
for (int j=0; j<undersat_gas_tables_[i][1].size(); ++j) {
undersat_gas_tables_[i][1][j] = 1.0/undersat_gas_tables_[i][1][j];
}
}
}

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

@ -471,8 +471,7 @@ namespace Opm
return cell_to_func_.empty() ? satfuncset_[0] : satfuncset_[cell_to_func_[cell]];
}
// Initialize saturation scaling parameter
// Initialize saturation scaling parameters
template <class SatFuncSet>
template <class T>
void SaturationPropsFromDeck<SatFuncSet>::initEPS(const EclipseGridParser& deck,
@ -1003,8 +1002,9 @@ namespace Opm
bool hasENPTVD = newParserDeck->hasKeyword("ENPTVD");
bool hasENKRVD = newParserDeck->hasKeyword("ENKRVD");
int itab = 0;
std::vector<std::vector<std::vector<double> > > table_dummy;
std::vector<std::vector<std::vector<double> > >& table = table_dummy;
std::vector<std::vector<double> > param_col;
std::vector<std::vector<double> > depth_col;
std::vector<std::string> col_names;
// Active keyword assigned default values for each cell (in case of possible box-wise assignment)
int phase_pos_aqua = phase_usage_.phase_pos[BlackoilPhases::Aqua];
@ -1070,22 +1070,14 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
Opm::EnptvdTable enptvd(newParserDeck->getKeyword("ENPTVD"));
table.resize(1); // only one region supported so far
for (unsigned i = 0; i < table.size(); ++i) {
table[i].resize(9);
for (int k = 0; k < enptvd.numRows(); ++k) {
table[i][0] = enptvd.getDepthColumn();
table[i][1] = enptvd.getSwcoColumn();
table[i][2] = enptvd.getSwcritColumn();
table[i][3] = enptvd.getSwmaxColumn();
table[i][4] = enptvd.getSgcoColumn();
table[i][5] = enptvd.getSgcritColumn();
table[i][6] = enptvd.getSgmaxColumn();
table[i][7] = enptvd.getSowcritColumn();
table[i][8] = enptvd.getSogcritColumn();
}
int num_tables = newParserDeck->getKeyword("ENPTVD")->size();
param_col.resize(num_tables);
depth_col.resize(num_tables);
col_names.resize(9);
for (int table_num=0; table_num<num_tables; ++table_num) {
Opm::SimpleTable enptvd(newParserDeck->getKeyword("ENPTVD"), col_names, table_num);
depth_col[table_num] = enptvd.getColumn(0); // depth
param_col[table_num] = enptvd.getColumn(itab); // itab=[1-8]: swl swcr swu sgl sgcr sgu sowcr sogcr
}
}
} else if ((keyword[0] == 'K' && (useKeyword || hasENKRVD)) || (keyword[1] == 'K' && useKeyword) ) {
@ -1141,21 +1133,14 @@ namespace Opm
OPM_THROW(std::runtime_error, " -- unknown keyword: '" << keyword << "'");
}
if (!useKeyword && itab > 0) {
Opm::EnkrvdTable enkrvd(newParserDeck->getKeyword("ENKRVD"));
table.resize(1); // only one region supported so far
for (unsigned i = 0; i < table.size(); ++i) {
table[i].resize(8);
for (int k = 0; k < enkrvd.numRows(); ++k) {
table[i][0] = enkrvd.getDepthColumn();
table[i][1] = enkrvd.getKrwmaxColumn();
table[i][2] = enkrvd.getKrgmaxColumn();
table[i][3] = enkrvd.getKromaxColumn();
table[i][4] = enkrvd.getKrwcritColumn();
table[i][5] = enkrvd.getKrgcritColumn();
table[i][6] = enkrvd.getKrocritgColumn();
table[i][7] = enkrvd.getKrocritwColumn();
}
int num_tables = newParserDeck->getKeyword("ENKRVD")->size();
param_col.resize(num_tables);
depth_col.resize(num_tables);
col_names.resize(8);
for (int table_num=0; table_num<num_tables; ++table_num) {
Opm::SimpleTable enkrvd(newParserDeck->getKeyword("ENKRVD"), col_names, table_num);
depth_col[table_num] = enkrvd.getColumn(0); // depth
param_col[table_num] = enkrvd.getColumn(itab); // itab=[1-7]: krw krg kro krwr krgr krorw krorg
}
}
}
@ -1172,25 +1157,15 @@ namespace Opm
scaleparam[c] = val[deck_pos];
}
} else {
// TODO for new parser
/*
std::cout << "--- Scaling parameter '" << keyword << "' assigned via ";
if (keyword[0] == 'S')
newParserDeck.getENPTVD().write(std::cout);
else
newParserDeck.getENKRVD().write(std::cout);
*/
const int dim = dimensions;
for (int cell = 0; cell < number_of_cells; ++cell) {
int jtab = cell_to_func_.empty() ? 0 : cell_to_func_[cell];
if (table[itab][jtab][0] != -1.0) {
std::vector<double>& depth = table[0][jtab];
std::vector<double>& val = table[itab][jtab];
if (param_col[jtab][0] >= 0.0) {
double zc = UgGridHelpers
::getCoordinate(UgGridHelpers::increment(begin_cell_centroid, cell, dim),
dim-1);
if (zc >= depth.front() && zc <= depth.back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth, val, zc);
if (zc >= depth_col[jtab].front() && zc <= depth_col[jtab].back()) { //don't want extrap outside depth interval
scaleparam[cell] = linearInterpolation(depth_col[jtab], param_col[jtab], zc);
}
}
}

829
opm/core/simulator/EquilibrationHelpers.hpp Normal file
View File

@ -0,0 +1,829 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#define OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/RootFinders.hpp>
#include <memory>
/*
---- synopsis of EquilibrationHelpers.hpp ----
namespace Opm
{
namespace Equil {
template <class Props>
class DensityCalculator;
template <>
class DensityCalculator< BlackoilPropertiesInterface >;
namespace Miscibility {
class RsFunction;
class NoMixing;
class RsVD;
class RsSatAtContact;
}
struct EquilRecord;
template <class DensCalc>
class EquilReg;
struct PcEq;
inline double satFromPc(const BlackoilPropertiesInterface& props,
const int phase,
const int cell,
const double target_pc,
const bool increasing = false);
struct PcEqSum
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
const int phase1,
const int phase2,
const int cell,
const double target_pc);
} // namespace Equil
} // namespace Opm
---- end of synopsis of EquilibrationHelpers.hpp ----
*/
namespace Opm
{
/**
* Types and routines that collectively implement a basic
* ECLIPSE-style equilibration-based initialisation scheme.
*
* This namespace is intentionally nested to avoid name clashes
* with other parts of OPM.
*/
namespace Equil {
template <class Props>
class DensityCalculator;
/**
* Facility for calculating phase densities based on the
* BlackoilPropertiesInterface.
*
* Implements the crucial <CODE>operator()(p,svol)</CODE>
* function that is expected by class EquilReg.
*/
template <>
class DensityCalculator< BlackoilPropertiesInterface > {
public:
/**
* Constructor.
*
* \param[in] props Implementation of the
* BlackoilPropertiesInterface.
*
* \param[in] c Single cell used as a representative cell
* in a PVT region.
*/
DensityCalculator(const BlackoilPropertiesInterface& props,
const int c)
: props_(props)
, c_(1, c)
{
}
/**
* Compute phase densities of all phases at phase point
* given by (pressure, surface volume) tuple.
*
* \param[in] p Fluid pressure.
*
* \param[in] z Surface volumes of all phases.
*
* \return Phase densities at phase point.
*/
std::vector<double>
operator()(const double p,
const std::vector<double>& z) const
{
const int np = props_.numPhases();
std::vector<double> A(np * np, 0);
assert (z.size() == std::vector<double>::size_type(np));
double* dAdp = 0;
props_.matrix(1, &p, &z[0], &c_[0], &A[0], dAdp);
std::vector<double> rho(np, 0.0);
props_.density(1, &A[0], &rho[0]);
return rho;
}
private:
const BlackoilPropertiesInterface& props_;
const std::vector<int> c_;
};
/**
* Types and routines relating to phase mixing in
* equilibration calculations.
*/
namespace Miscibility {
/**
* Base class for phase mixing functions.
*/
class RsFunction
{
public:
/**
* Function call operator.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
virtual double operator()(const double depth,
const double press,
const double sat = 0.0) const = 0;
};
/**
* Type that implements "no phase mixing" policy.
*/
class NoMixing : public RsFunction {
public:
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press. In "no mixing
* policy", this is identically zero.
*/
double
operator()(const double /* depth */,
const double /* press */,
const double sat = 0.0) const
{
return 0.0;
}
};
/**
* Type that implements "dissolved gas-oil ratio"
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RSVD'.
*/
class RsVD : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] depth Depth nodes.
* \param[in] rs Dissolved gas-oil ratio at @c depth.
*/
RsVD(const BlackoilPropertiesInterface& props,
const int cell,
const std::vector<double>& depth,
const std::vector<double>& rs)
: props_(props)
, cell_(cell)
, depth_(depth)
, rs_(rs)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double depth,
const double press,
const double sat_gas = 0.0) const
{
if (sat_gas > 0.0) {
return satRs(press);
} else {
return std::min(satRs(press), linearInterpolation(depth_, rs_, depth));
}
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
std::vector<double> depth_; /**< Depth nodes */
std::vector<double> rs_; /**< Dissolved gas-oil ratio */
double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*opos + gpos] / A_[np*opos + opos];
}
};
/**
* Type that implements "vaporized oil-gas ratio"
* tabulated as a function of depth policy. Data
* typically taken from keyword 'RVVD'.
*/
class RvVD : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] depth Depth nodes.
* \param[in] rv Dissolved gas-oil ratio at @c depth.
*/
RvVD(const BlackoilPropertiesInterface& props,
const int cell,
const std::vector<double>& depth,
const std::vector<double>& rv)
: props_(props)
, cell_(cell)
, depth_(depth)
, rv_(rv)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RV
* value.
*
* \param[in] press Pressure at which to calculate RV
* value.
*
* \return Vaporized oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double depth,
const double press,
const double sat_oil = 0.0 ) const
{
if (sat_oil > 0.0) {
return satRv(press);
} else {
return std::min(satRv(press), linearInterpolation(depth_, rv_, depth));
}
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
std::vector<double> depth_; /**< Depth nodes */
std::vector<double> rv_; /**< Vaporized oil-gas ratio */
double z_[BlackoilPhases::MaxNumPhases];
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*gpos + opos] / A_[np*gpos + gpos];
}
};
/**
* Class that implements "dissolved gas-oil ratio" (Rs)
* as function of depth and pressure as follows:
*
* 1. The Rs at the gas-oil contact is equal to the
* saturated Rs value, Rs_sat_contact.
*
* 2. The Rs elsewhere is equal to Rs_sat_contact, but
* constrained to the saturated value as given by the
* local pressure.
*
* This should yield Rs-values that are constant below the
* contact, and decreasing above the contact.
*/
class RsSatAtContact : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] p_contact oil pressure at the contact
*/
RsSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
: props_(props), cell_(cell)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1e100;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0;
rs_sat_contact_ = satRs(p_contact);
}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RS
* value.
*
* \param[in] press Pressure at which to calculate RS
* value.
*
* \return Dissolved gas-oil ratio (RS) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double /* depth */,
const double press,
const double sat_gas = 0.0) const
{
if (sat_gas > 0.0) {
return satRs(press);
} else {
return std::min(satRs(press), rs_sat_contact_);
}
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
double z_[BlackoilPhases::MaxNumPhases];
double rs_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRs(const double press) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
// Rs/Bo is in the gas row and oil column of A_.
// 1/Bo is in the oil row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*opos + gpos] / A_[np*opos + opos];
}
};
/**
* Class that implements "vaporized oil-gas ratio" (Rv)
* as function of depth and pressure as follows:
*
* 1. The Rv at the gas-oil contact is equal to the
* saturated Rv value, Rv_sat_contact.
*
* 2. The Rv elsewhere is equal to Rv_sat_contact, but
* constrained to the saturated value as given by the
* local pressure.
*
* This should yield Rv-values that are constant below the
* contact, and decreasing above the contact.
*/
class RvSatAtContact : public RsFunction {
public:
/**
* Constructor.
*
* \param[in] props property object
* \param[in] cell any cell in the pvt region
* \param[in] p_contact oil pressure at the contact
*/
RvSatAtContact(const BlackoilPropertiesInterface& props, const int cell, const double p_contact)
: props_(props), cell_(cell)
{
auto pu = props_.phaseUsage();
std::fill(z_, z_ + BlackoilPhases::MaxNumPhases, 0.0);
z_[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0;
z_[pu.phase_pos[BlackoilPhases::Liquid]] = 1e100;
rv_sat_contact_ = satRv(p_contact);
}
/**
* Function call.
*
* \param[in] depth Depth at which to calculate RV
* value.
*
* \param[in] press Pressure at which to calculate RV
* value.
*
* \return Dissolved oil-gas ratio (RV) at depth @c
* depth and pressure @c press.
*/
double
operator()(const double /*depth*/,
const double press,
const double sat_oil = 0.0) const
{
if (sat_oil > 0.0) {
return satRv(press);
} else {
return std::min(satRv(press), rv_sat_contact_);
}
}
private:
const BlackoilPropertiesInterface& props_;
const int cell_;
double z_[BlackoilPhases::MaxNumPhases];
double rv_sat_contact_;
mutable double A_[BlackoilPhases::MaxNumPhases * BlackoilPhases::MaxNumPhases];
double satRv(const double press) const
{
props_.matrix(1, &press, z_, &cell_, A_, 0);
// Rv/Bg is in the oil row and gas column of A_.
// 1/Bg is in the gas row and column.
// Recall also that it is stored in column-major order.
const int opos = props_.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gpos = props_.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const int np = props_.numPhases();
return A_[np*gpos + opos] / A_[np*gpos + gpos];
}
};
} // namespace Miscibility
/**
* Equilibration record.
*
* Layout and contents inspired by first six items of
* ECLIPSE's 'EQUIL' records. This is the minimum amount of
* input data needed to define phase pressures in an
* equilibration region.
*
* Data consists of three pairs of depth and pressure values:
* 1. main
* - @c depth Main datum depth.
* - @c press Pressure at datum depth.
*
* 2. woc
* - @c depth Depth of water-oil contact
* - @c press water-oil capillary pressure at water-oil contact.
* Capillary pressure defined as "P_oil - P_water".
*
* 3. goc
* - @c depth Depth of gas-oil contact
* - @c press Gas-oil capillary pressure at gas-oil contact.
* Capillary pressure defined as "P_gas - P_oil".
*/
struct EquilRecord {
struct {
double depth;
double press;
} main, woc, goc;
int live_oil_table_index;
int wet_gas_table_index;
int N;
};
/**
* Aggregate information base of an equilibration region.
*
* Provides inquiry methods for retrieving depths of contacs
* and pressure values as well as a means of calculating fluid
* densities, dissolved gas-oil ratio and vapourised oil-gas
* ratios.
*
* \tparam DensCalc Type that provides access to a phase
* density calculation facility. Must implement an operator()
* declared as
* <CODE>
* std::vector<double>
* operator()(const double press,
* const std::vector<double>& svol )
* </CODE>
* that calculates the phase densities of all phases in @c
* svol at fluid pressure @c press.
*/
template <class DensCalc>
class EquilReg {
public:
/**
* Constructor.
*
* \param[in] rec Equilibration data of current region.
* \param[in] density Density calculator of current region.
* \param[in] rs Calculator of dissolved gas-oil ratio.
* \param[in] rv Calculator of vapourised oil-gas ratio.
* \param[in] pu Summary of current active phases.
*/
EquilReg(const EquilRecord& rec,
const DensCalc& density,
std::shared_ptr<Miscibility::RsFunction> rs,
std::shared_ptr<Miscibility::RsFunction> rv,
const PhaseUsage& pu)
: rec_ (rec)
, density_(density)
, rs_ (rs)
, rv_ (rv)
, pu_ (pu)
{
}
/**
* Type of density calculator.
*/
typedef DensCalc CalcDensity;
/**
* Type of dissolved gas-oil ratio calculator.
*/
typedef Miscibility::RsFunction CalcDissolution;
/**
* Type of vapourised oil-gas ratio calculator.
*/
typedef Miscibility::RsFunction CalcEvaporation;
/**
* Datum depth in current region
*/
double datum() const { return this->rec_.main.depth; }
/**
* Pressure at datum depth in current region.
*/
double pressure() const { return this->rec_.main.press; }
/**
* Depth of water-oil contact.
*/
double zwoc() const { return this->rec_.woc .depth; }
/**
* water-oil capillary pressure at water-oil contact.
*
* \return P_o - P_w at WOC.
*/
double pcow_woc() const { return this->rec_.woc .press; }
/**
* Depth of gas-oil contact.
*/
double zgoc() const { return this->rec_.goc .depth; }
/**
* Gas-oil capillary pressure at gas-oil contact.
*
* \return P_g - P_o at GOC.
*/
double pcgo_goc() const { return this->rec_.goc .press; }
/**
* Retrieve phase density calculator of current region.
*/
const CalcDensity&
densityCalculator() const { return this->density_; }
/**
* Retrieve dissolved gas-oil ratio calculator of current
* region.
*/
const CalcDissolution&
dissolutionCalculator() const { return *this->rs_; }
/**
* Retrieve vapourised oil-gas ratio calculator of current
* region.
*/
const CalcEvaporation&
evaporationCalculator() const { return *this->rv_; }
/**
* Retrieve active fluid phase summary.
*/
const PhaseUsage&
phaseUsage() const { return this->pu_; }
private:
EquilRecord rec_; /**< Equilibration data */
DensCalc density_; /**< Density calculator */
std::shared_ptr<Miscibility::RsFunction> rs_; /**< RS calculator */
std::shared_ptr<Miscibility::RsFunction> rv_; /**< RV calculator */
PhaseUsage pu_; /**< Active phase summary */
};
/// Functor for inverting capillary pressure function.
/// Function represented is
/// f(s) = pc(s) - target_pc
struct PcEq
{
PcEq(const BlackoilPropertiesInterface& props,
const int phase,
const int cell,
const double target_pc)
: props_(props),
phase_(phase),
cell_(cell),
target_pc_(target_pc)
{
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
}
double operator()(double s) const
{
s_[phase_] = s;
props_.capPress(1, s_, &cell_, pc_, 0);
return pc_[phase_] - target_pc_;
}
private:
const BlackoilPropertiesInterface& props_;
const int phase_;
const int cell_;
const int target_pc_;
mutable double s_[BlackoilPhases::MaxNumPhases];
mutable double pc_[BlackoilPhases::MaxNumPhases];
};
/// Compute saturation of some phase corresponding to a given
/// capillary pressure.
inline double satFromPc(const BlackoilPropertiesInterface& props,
const int phase,
const int cell,
const double target_pc,
const bool increasing = false)
{
// Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases];
double smaxarr[BlackoilPhases::MaxNumPhases];
props.satRange(1, &cell, sminarr, smaxarr);
const double s0 = increasing ? smaxarr[phase] : sminarr[phase];
const double s1 = increasing ? sminarr[phase] : smaxarr[phase];
// Create the equation f(s) = pc(s) - target_pc
const PcEq f(props, phase, cell, target_pc);
const double f0 = f(s0);
const double f1 = f(s1);
if (f0 <= 0.0) {
return s0;
} else if (f1 > 0.0) {
return s1;
} else {
const int max_iter = 30;
const double tol = 1e-6;
int iter_used = -1;
typedef RegulaFalsi<ThrowOnError> ScalarSolver;
const double sol = ScalarSolver::solve(f, std::min(s0, s1), std::max(s0, s1), max_iter, tol, iter_used);
return sol;
}
}
/// Functor for inverting a sum of capillary pressure functions.
/// Function represented is
/// f(s) = pc1(s) + pc2(1 - s) - target_pc
struct PcEqSum
{
PcEqSum(const BlackoilPropertiesInterface& props,
const int phase1,
const int phase2,
const int cell,
const double target_pc)
: props_(props),
phase1_(phase1),
phase2_(phase2),
cell_(cell),
target_pc_(target_pc)
{
std::fill(s_, s_ + BlackoilPhases::MaxNumPhases, 0.0);
std::fill(pc_, pc_ + BlackoilPhases::MaxNumPhases, 0.0);
}
double operator()(double s) const
{
s_[phase1_] = s;
s_[phase2_] = 1.0 - s;
props_.capPress(1, s_, &cell_, pc_, 0);
return pc_[phase1_] + pc_[phase2_] - target_pc_;
}
private:
const BlackoilPropertiesInterface& props_;
const int phase1_;
const int phase2_;
const int cell_;
const int target_pc_;
mutable double s_[BlackoilPhases::MaxNumPhases];
mutable double pc_[BlackoilPhases::MaxNumPhases];
};
/// Compute saturation of some phase corresponding to a given
/// capillary pressure, where the capillary pressure function
/// is given as a sum of two other functions.
inline double satFromSumOfPcs(const BlackoilPropertiesInterface& props,
const int phase1,
const int phase2,
const int cell,
const double target_pc)
{
// Find minimum and maximum saturations.
double sminarr[BlackoilPhases::MaxNumPhases];
double smaxarr[BlackoilPhases::MaxNumPhases];
props.satRange(1, &cell, sminarr, smaxarr);
const double smin = sminarr[phase1];
const double smax = smaxarr[phase1];
// Create the equation f(s) = pc1(s) + pc2(1-s) - target_pc
const PcEqSum f(props, phase1, phase2, cell, target_pc);
const double f0 = f(smin);
const double f1 = f(smax);
if (f0 <= 0.0) {
return smin;
} else if (f1 > 0.0) {
return smax;
} else {
const int max_iter = 30;
const double tol = 1e-6;
int iter_used = -1;
typedef RegulaFalsi<ThrowOnError> ScalarSolver;
const double sol = ScalarSolver::solve(f, smin, smax, max_iter, tol, iter_used);
return sol;
}
}
} // namespace Equil
} // namespace Opm
#endif // OPM_EQUILIBRATIONHELPERS_HEADER_INCLUDED

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

@ -20,7 +20,8 @@
#include "SimulatorOutput.hpp"
// we need complete definitions for these types
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/TimeMap.hpp>
#include <opm/core/io/OutputWriter.hpp>
#include <opm/core/simulator/SimulatorTimer.hpp>
@ -30,7 +31,8 @@ using namespace Opm;
SimulatorOutputBase::SimulatorOutputBase (
const parameter::ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const Deck> parser,
std::shared_ptr <const TimeMap> timeMap,
std::shared_ptr <const UnstructuredGrid> grid,
std::shared_ptr <const SimulatorTimer> timer,
std::shared_ptr <const SimulatorState> state,
@ -39,6 +41,7 @@ SimulatorOutputBase::SimulatorOutputBase (
// store all parameters passed into the object, making them curried
// parameters to the writeOutput function.
: timer_ (timer )
, timeMap_ (timeMap )
, reservoirState_ (state )
, wellState_ (wellState)
@ -49,16 +52,8 @@ SimulatorOutputBase::SimulatorOutputBase (
// always start from the first timestep
, next_ (0) {
// make a list of times to dump. since the original list are relative
// timesteps, we make a list of accumulated such to compare with
// current time. add an extra zero at the beginning so that the
// initial state is also written
const std::vector <double>& tstep = parser->getTSTEP ().tstep_;
times_.resize (tstep.size (), 0.);
std::partial_sum (tstep.begin(), tstep.end(), times_.begin());
// write the static initialization files, even before simulation starts
writer_->writeInit (*timer, *state, *wellState);
writer_->writeInit (*timer);
}
// default destructor is OK, just need to be defined
@ -76,15 +71,18 @@ SimulatorOutputBase::writeOutput () {
// if the simulator signals for timesteps that aren't reporting
// times, then ignore them
if (next_ < times_.size () && times_[next_] <= this_time) {
if (next_ < timeMap_->size ()
&& timeMap_->getTimePassedUntil (next_) <= this_time) {
// uh-oh, the simulator has skipped reporting timesteps that
// occurred before this timestep (it doesn't honor the TSTEP setting)
while (next_ < times_.size () && times_[next_] < this_time) {
while (next_ < timeMap_->size ()
&& timeMap_->getTimePassedUntil (next_) < this_time) {
++next_;
}
// report this timestep if it matches
if (next_ < times_.size () && times_[next_] == this_time) {
if (next_ < timeMap_->size ()
&& timeMap_->getTimePassedUntil (next_) == this_time) {
// make sure the simulator has spilled all necessary internal
// state. notice that this calls *our* sync, which is overridden
// in the template companion to call the simulator

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

@ -32,11 +32,12 @@ struct UnstructuredGrid;
namespace Opm {
// forward definitions
class EclipseGridParser;
class Deck;
class OutputWriter;
namespace parameter { class ParameterGroup; }
class SimulatorState;
class SimulatorTimer;
class TimeMap;
class WellState;
/**
@ -53,7 +54,8 @@ protected:
* need to pick them up from the object members.
*/
SimulatorOutputBase (const parameter::ParameterGroup& p,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const Deck> parser,
std::shared_ptr <const TimeMap> timeMap,
std::shared_ptr <const UnstructuredGrid> grid,
std::shared_ptr <const SimulatorTimer> timer,
std::shared_ptr <const SimulatorState> state,
@ -75,6 +77,7 @@ protected:
/// Just hold a reference to these objects that are owned elsewhere.
std::shared_ptr <const SimulatorTimer> timer_;
std::shared_ptr <const TimeMap> timeMap_;
std::shared_ptr <const SimulatorState> reservoirState_;
std::shared_ptr <const WellState> wellState_;
@ -110,7 +113,7 @@ private:
* ParameterGroup params (argc, argv, false);
*
* // input file
* auto deck = make_shared <EclipseGridParser> ( ... );
* auto deck = make_shared <const Deck> ( ... );
* const GridManager manager (*parser);
* auto grid = share_obj (*manager.c_grid ());
*
@ -122,11 +125,12 @@ private:
* auto wellState = make_shared <WellState> ();
*
* // set up simulation
* auto timeMap = make_shared <const TimeMap> (deck);
* auto sim = make_shared <SimulatorIncompTwophase> (params, *grid, ... );
*
* // use this to dump state to disk
* auto output = make_shared <SimulatorOutput> (
* params, deck, grid, timer, state, wellState, sim);
* params, deck, timeMap, grid, timer, state, wellState, sim);
*
* // start simulation
* sim.run (timer, state, ... )
@ -139,14 +143,16 @@ private:
template <typename Simulator>
struct SimulatorOutput : public SimulatorOutputBase {
SimulatorOutput (const parameter::ParameterGroup& params,
std::shared_ptr <EclipseGridParser> parser,
std::shared_ptr <const Deck> parser,
std::shared_ptr <const TimeMap> timeMap,
std::shared_ptr <const UnstructuredGrid> grid,
std::shared_ptr <const SimulatorTimer> timer,
std::shared_ptr <const SimulatorState> state,
std::shared_ptr <const WellState> wellState,
std::shared_ptr <Simulator> sim)
// send all other parameters to base class
: SimulatorOutputBase (params, parser, grid, timer, state, wellState)
: SimulatorOutputBase (params, parser, timeMap,
grid, timer, state, wellState)
// store reference to simulator in derived class
, sim_ (sim) {
@ -161,7 +167,8 @@ struct SimulatorOutput : public SimulatorOutputBase {
* the arguments passed exceeds the lifetime of this object.
*/
SimulatorOutput (const parameter::ParameterGroup& params,
EclipseGridParser& parser,
const Deck& parser,
const TimeMap& timeMap,
const UnstructuredGrid& grid,
const SimulatorTimer& timer,
const SimulatorState& state,
@ -170,6 +177,7 @@ struct SimulatorOutput : public SimulatorOutputBase {
// send all other parameters to base class
: SimulatorOutputBase (params,
share_obj (parser),
share_obj (timeMap),
share_obj (grid),
share_obj (timer),
share_obj (state),

673
opm/core/simulator/initStateEquil.hpp Normal file
View File

@ -0,0 +1,673 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_INITSTATEEQUIL_HEADER_INCLUDED
#define OPM_INITSTATEEQUIL_HEADER_INCLUDED
#include <opm/core/simulator/EquilibrationHelpers.hpp>
#include <opm/core/simulator/BlackoilState.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPropertiesInterface.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/core/utility/RegionMapping.hpp>
#include <opm/core/utility/Units.hpp>
#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
#include <opm/parser/eclipse/Utility/SimpleTable.hpp>
#include <array>
#include <cassert>
#include <utility>
#include <vector>
/**
* \file
* Facilities for an ECLIPSE-style equilibration-based
* initialisation scheme (keyword 'EQUIL').
*/
struct UnstructuredGrid;
namespace Opm
{
/**
* Compute initial state by an equilibration procedure.
*
* The following state fields are modified:
* pressure(),
* saturation(),
* surfacevol(),
* gasoilratio(),
* rv().
*
* \param[in] grid Grid.
* \param[in] props Property object, pvt and capillary properties are used.
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck,
const double gravity,
BlackoilState& state);
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const double gravity,
BlackoilState& state);
/**
* Types and routines that collectively implement a basic
* ECLIPSE-style equilibration-based initialisation scheme.
*
* This namespace is intentionally nested to avoid name clashes
* with other parts of OPM.
*/
namespace Equil {
/**
* Compute initial phase pressures by means of equilibration.
*
* This function uses the information contained in an
* equilibration record (i.e., depths and pressurs) as well as
* a density calculator and related data to vertically
* integrate the phase pressure ODE
* \f[
* \frac{\mathrm{d}p_{\alpha}}{\mathrm{d}z} =
* \rho_{\alpha}(z,p_{\alpha})\cdot g
* \f]
* in which \f$\rho_{\alpha}$ denotes the fluid density of
* fluid phase \f$\alpha\f$, \f$p_{\alpha}\f$ is the
* corresponding phase pressure, \f$z\f$ is the depth and
* \f$g\f$ is the acceleration due to gravity (assumed
* directed downwords, in the positive \f$z\f$ direction).
*
* \tparam Region Type of an equilibration region information
* base. Typically an instance of the EquilReg
* class template.
*
* \tparam CellRange Type of cell range that demarcates the
* cells pertaining to the current
* equilibration region. Must implement
* methods begin() and end() to bound the range
* as well as provide an inner type,
* const_iterator, to traverse the range.
*
* \param[in] G Grid.
* \param[in] reg Current equilibration region.
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] grav Acceleration of gravity.
*
* \return Phase pressures, one vector for each active phase,
* of pressure values in each cell in the current
* equilibration region.
*/
template <class Region, class CellRange>
std::vector< std::vector<double> >
phasePressures(const UnstructuredGrid& G,
const Region& reg,
const CellRange& cells,
const double grav = unit::gravity);
/**
* Compute initial phase saturations by means of equilibration.
*
* \tparam Region Type of an equilibration region information
* base. Typically an instance of the EquilReg
* class template.
*
* \tparam CellRange Type of cell range that demarcates the
* cells pertaining to the current
* equilibration region. Must implement
* methods begin() and end() to bound the range
* as well as provide an inner type,
* const_iterator, to traverse the range.
*
* \param[in] reg Current equilibration region.
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] props Property object, needed for capillary functions.
* \param[in] phase_pressures Phase pressures, one vector for each active phase,
* of pressure values in each cell in the current
* equilibration region.
* \return Phase saturations, one vector for each phase, each containing
* one saturation value per cell in the region.
*/
template <class Region, class CellRange>
std::vector< std::vector<double> >
phaseSaturations(const Region& reg,
const CellRange& cells,
const BlackoilPropertiesInterface& props,
const std::vector< std::vector<double> >& phase_pressures);
/**
* Compute initial Rs values.
*
* \tparam CellRangeType Type of cell range that demarcates the
* cells pertaining to the current
* equilibration region. Must implement
* methods begin() and end() to bound the range
* as well as provide an inner type,
* const_iterator, to traverse the range.
*
* \param[in] grid Grid.
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] oil_pressure Oil pressure for each cell in range.
* \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range.
*/
template <class CellRangeType>
std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells,
const std::vector<double> oil_pressure,
const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation);
namespace DeckDependent {
inline
std::vector<EquilRecord>
getEquil(const EclipseGridParser& deck)
{
if (deck.hasField("EQUIL")) {
const EQUIL& eql = deck.getEQUIL();
typedef std::vector<EquilLine>::size_type sz_t;
const sz_t nrec = eql.equil.size();
std::vector<EquilRecord> ret;
ret.reserve(nrec);
for (sz_t r = 0; r < nrec; ++r) {
const EquilLine& rec = eql.equil[r];
EquilRecord record =
{
{ rec.datum_depth_ ,
rec.datum_depth_pressure_ }
,
{ rec.water_oil_contact_depth_ ,
rec.oil_water_cap_pressure_ }
,
{ rec.gas_oil_contact_depth_ ,
rec.gas_oil_cap_pressure_ }
,
rec.live_oil_table_index_
,
rec.wet_gas_table_index_
,
rec.N_
};
if (record.N != 0) {
OPM_THROW(std::domain_error,
"kw EQUIL, item 9: Only N=0 supported.");
}
ret.push_back(record);
}
return ret;
}
else {
OPM_THROW(std::domain_error,
"Deck does not provide equilibration data.");
}
}
inline
std::vector<EquilRecord>
getEquil(const Opm::DeckConstPtr newParserDeck)
{
if (newParserDeck->hasKeyword("EQUIL")) {
Opm::EquilWrapper eql(newParserDeck->getKeyword("EQUIL"));
const int nrec = eql.numRegions();
std::vector<EquilRecord> ret;
ret.reserve(nrec);
for (int r = 0; r < nrec; ++r) {
EquilRecord record =
{
{ eql.datumDepth(r) ,
eql.datumDepthPressure(r) }
,
{ eql.waterOilContactDepth(r) ,
eql.waterOilContactCapillaryPressure(r) }
,
{ eql.gasOilContactDepth(r) ,
eql.gasOilContactCapillaryPressure(r) }
,
eql.liveOilInitProceedure(r)
,
eql.wetGasInitProceedure(r)
,
eql.initializationTargetAccuracy(r)
};
if (record.N != 0) {
OPM_THROW(std::domain_error,
"kw EQUIL, item 9: Only N=0 supported.");
}
ret.push_back(record);
}
return ret;
}
else {
OPM_THROW(std::domain_error,
"Deck does not provide equilibration data.");
}
}
inline
std::vector<int>
equilnum(const EclipseGridParser& deck,
const UnstructuredGrid& G )
{
std::vector<int> eqlnum;
if (deck.hasField("EQLNUM")) {
const std::vector<int>& e = deck.getIntegerValue("EQLNUM");
eqlnum.reserve(e.size());
std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
std::bind2nd(std::minus<int>(), 1));
}
else {
// No explicit equilibration region.
// All cells in region zero.
eqlnum.assign(G.number_of_cells, 0);
}
return eqlnum;
}
inline
std::vector<int>
equilnum(const Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& G )
{
std::vector<int> eqlnum;
if (newParserDeck->hasKeyword("EQLNUM")) {
const std::vector<int>& e = newParserDeck->getKeyword("EQLNUM")->getIntData();
eqlnum.reserve(e.size());
std::transform(e.begin(), e.end(), std::back_inserter(eqlnum),
std::bind2nd(std::minus<int>(), 1));
}
else {
// No explicit equilibration region.
// All cells in region zero.
eqlnum.assign(G.number_of_cells, 0);
}
return eqlnum;
}
template <class InputDeck>
class InitialStateComputer;
template <>
class InitialStateComputer<Opm::EclipseGridParser> {
public:
InitialStateComputer(const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck ,
const UnstructuredGrid& G ,
const double grav = unit::gravity)
: pp_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
sat_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
rs_(G.number_of_cells),
rv_(G.number_of_cells)
{
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(deck);
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(deck, G));
// Create Rs functions.
rs_func_.reserve(rec.size());
if (deck.hasField("DISGAS")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].live_oil_table_index > 0) {
if (deck.hasField("RSVD")) {
// TODO When this kw is actually parsed, also check for proper number of available tables
// For now, just use dummy ...
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
std::vector<double> rs; rs.push_back(0.0); rs.push_back(100.0);
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, depth, rs));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RSVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press;
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
rv_func_.reserve(rec.size());
if (deck.hasField("VAPOIL")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].wet_gas_table_index > 0) {
if (deck.hasField("RVVD")) {
// TODO When this kw is actually parsed, also check for proper number of available tables
// For now, just use dummy ...
std::vector<double> depth; depth.push_back(0.0); depth.push_back(100.0);
std::vector<double> rv; rv.push_back(0.0); rv.push_back(0.0001);
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, depth, rv));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RVVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press + rec[i].goc.press;
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, props, G, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
}
typedef std::vector<double> Vec;
typedef std::vector<Vec> PVec; // One per phase.
const PVec& press() const { return pp_; }
const PVec& saturation() const { return sat_; }
const Vec& rs() const { return rs_; }
const Vec& rv() const { return rv_; }
private:
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
typedef EquilReg<RhoCalc> EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
PVec pp_;
PVec sat_;
Vec rs_;
Vec rv_;
template <class RMap>
void
calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec ,
const Opm::BlackoilPropertiesInterface& props,
const UnstructuredGrid& G ,
const double grav)
{
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)
{
const typename RMap::CellRange cells = reg.cells(r);
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell);
const EqReg eqreg(rec[r], calc,
rs_func_[r], rv_func_[r],
props.phaseUsage());
PVec press = phasePressures(G, eqreg, cells, grav);
const PVec sat = phaseSaturations(eqreg, cells, props, press);
const int np = props.numPhases();
for (int p = 0; p < np; ++p) {
copyFromRegion(press[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]);
}
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_);
}
}
}
template <class CellRangeType>
void copyFromRegion(const Vec& source,
const CellRangeType& cells,
Vec& destination)
{
auto s = source.begin();
auto c = cells.begin();
const auto e = cells.end();
for (; c != e; ++c, ++s) {
destination[*c] = *s;
}
}
};
template <>
class InitialStateComputer<Opm::DeckConstPtr> {
public:
InitialStateComputer(const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& G ,
const double grav = unit::gravity)
: pp_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
sat_(props.numPhases(),
std::vector<double>(G.number_of_cells)),
rs_(G.number_of_cells),
rv_(G.number_of_cells)
{
// Get the equilibration records.
const std::vector<EquilRecord> rec = getEquil(newParserDeck);
// Create (inverse) region mapping.
const RegionMapping<> eqlmap(equilnum(newParserDeck, G));
// Create Rs functions.
rs_func_.reserve(rec.size());
if (newParserDeck->hasKeyword("DISGAS")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].live_oil_table_index > 0) {
if (newParserDeck->hasKeyword("RSVD") && rec[i].live_oil_table_index <= newParserDeck->getKeyword("RSVD")->size()) {
Opm::SimpleTable rsvd(newParserDeck->getKeyword("RSVD"),std::vector<std::string>{"vd", "rs"},rec[i].live_oil_table_index-1);
std::vector<double> vd(rsvd.getColumn("vd"));
std::vector<double> rs(rsvd.getColumn("rs"));
rs_func_.push_back(std::make_shared<Miscibility::RsVD>(props, cell, vd, rs));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RSVD table " << (rec[i].live_oil_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RSVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press;
rs_func_.push_back(std::make_shared<Miscibility::RsSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rs_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
rv_func_.reserve(rec.size());
if (newParserDeck->hasKeyword("VAPOIL")) {
for (size_t i = 0; i < rec.size(); ++i) {
const int cell = *(eqlmap.cells(i).begin());
if (rec[i].wet_gas_table_index > 0) {
if (newParserDeck->hasKeyword("RVVD") && rec[i].wet_gas_table_index <= newParserDeck->getKeyword("RVVD")->size()) {
Opm::SimpleTable rvvd(newParserDeck->getKeyword("RVVD"),std::vector<std::string>{"vd", "rv"},rec[i].wet_gas_table_index-1);
std::vector<double> vd(rvvd.getColumn("vd"));
std::vector<double> rv(rvvd.getColumn("rv"));
rv_func_.push_back(std::make_shared<Miscibility::RvVD>(props, cell, vd, rv));
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: RVVD table " << (rec[i].wet_gas_table_index) << " not available.");
}
} else {
if (rec[i].goc.depth != rec[i].main.depth) {
OPM_THROW(std::runtime_error,
"Cannot initialise: when no explicit RVVD table is given, \n"
"datum depth must be at the gas-oil-contact. "
"In EQUIL region " << (i + 1) << " (counting from 1), this does not hold.");
}
const double p_contact = rec[i].main.press + rec[i].goc.press;
rv_func_.push_back(std::make_shared<Miscibility::RvSatAtContact>(props, cell, p_contact));
}
}
} else {
for (size_t i = 0; i < rec.size(); ++i) {
rv_func_.push_back(std::make_shared<Miscibility::NoMixing>());
}
}
// Compute pressures, saturations, rs and rv factors.
calcPressSatRsRv(eqlmap, rec, props, G, grav);
// Modify oil pressure in no-oil regions so that the pressures of present phases can
// be recovered from the oil pressure and capillary relations.
}
typedef std::vector<double> Vec;
typedef std::vector<Vec> PVec; // One per phase.
const PVec& press() const { return pp_; }
const PVec& saturation() const { return sat_; }
const Vec& rs() const { return rs_; }
const Vec& rv() const { return rv_; }
private:
typedef DensityCalculator<BlackoilPropertiesInterface> RhoCalc;
typedef EquilReg<RhoCalc> EqReg;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rs_func_;
std::vector< std::shared_ptr<Miscibility::RsFunction> > rv_func_;
PVec pp_;
PVec sat_;
Vec rs_;
Vec rv_;
template <class RMap>
void
calcPressSatRsRv(const RMap& reg ,
const std::vector< EquilRecord >& rec ,
const Opm::BlackoilPropertiesInterface& props,
const UnstructuredGrid& G ,
const double grav)
{
typedef Miscibility::NoMixing NoMix;
for (typename RMap::RegionId
r = 0, nr = reg.numRegions();
r < nr; ++r)
{
const typename RMap::CellRange cells = reg.cells(r);
const int repcell = *cells.begin();
const RhoCalc calc(props, repcell);
const EqReg eqreg(rec[r], calc,
rs_func_[r], rv_func_[r],
props.phaseUsage());
PVec press = phasePressures(G, eqreg, cells, grav);
const PVec sat = phaseSaturations(eqreg, cells, props, press);
const int np = props.numPhases();
for (int p = 0; p < np; ++p) {
copyFromRegion(press[p], cells, pp_[p]);
copyFromRegion(sat[p], cells, sat_[p]);
}
if (props.phaseUsage().phase_used[BlackoilPhases::Liquid]
&& props.phaseUsage().phase_used[BlackoilPhases::Vapour]) {
const int oilpos = props.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int gaspos = props.phaseUsage().phase_pos[BlackoilPhases::Vapour];
const Vec rs = computeRs(G, cells, press[oilpos], *(rs_func_[r]), sat[gaspos]);
const Vec rv = computeRs(G, cells, press[gaspos], *(rv_func_[r]), sat[oilpos]);
copyFromRegion(rs, cells, rs_);
copyFromRegion(rv, cells, rv_);
}
}
}
template <class CellRangeType>
void copyFromRegion(const Vec& source,
const CellRangeType& cells,
Vec& destination)
{
auto s = source.begin();
auto c = cells.begin();
const auto e = cells.end();
for (; c != e; ++c, ++s) {
destination[*c] = *s;
}
}
};
} // namespace DeckDependent
} // namespace Equil
} // namespace Opm
#include <opm/core/simulator/initStateEquil_impl.hpp>
#endif // OPM_INITSTATEEQUIL_HEADER_INCLUDED

781
opm/core/simulator/initStateEquil_impl.hpp Normal file
View File

@ -0,0 +1,781 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
#define OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED
#include <opm/core/grid.h>
#include <opm/core/props/BlackoilPhases.hpp>
#include <cassert>
#include <cmath>
#include <functional>
#include <vector>
namespace Opm
{
namespace Details {
template <class RHS>
class RK4IVP {
public:
RK4IVP(const RHS& f ,
const std::array<double,2>& span,
const double y0 ,
const int N )
: N_(N)
, span_(span)
{
const double h = stepsize();
const double h2 = h / 2;
const double h6 = h / 6;
y_.reserve(N + 1);
f_.reserve(N + 1);
y_.push_back(y0);
f_.push_back(f(span_[0], y0));
for (int i = 0; i < N; ++i) {
const double x = span_[0] + i*h;
const double y = y_.back();
const double k1 = f_[i];
const double k2 = f(x + h2, y + h2*k1);
const double k3 = f(x + h2, y + h2*k2);
const double k4 = f(x + h , y + h*k3);
y_.push_back(y + h6*(k1 + 2*(k2 + k3) + k4));
f_.push_back(f(x + h, y_.back()));
}
assert (y_.size() == std::vector<double>::size_type(N + 1));
}
double
operator()(const double x) const
{
// Dense output (O(h**3)) according to Shampine
// (Hermite interpolation)
const double h = stepsize();
int i = (x - span_[0]) / h;
const double t = (x - (span_[0] + i*h)) / h;
// Crude handling of evaluation point outside "span_";
if (i < 0) { i = 0; }
if (N_ <= i) { i = N_ - 1; }
const double y0 = y_[i], y1 = y_[i + 1];
const double f0 = f_[i], f1 = f_[i + 1];
double u = (1 - 2*t) * (y1 - y0);
u += h * ((t - 1)*f0 + t*f1);
u *= t * (t - 1);
u += (1 - t)*y0 + t*y1;
return u;
}
private:
int N_;
std::array<double,2> span_;
std::vector<double> y_;
std::vector<double> f_;
double
stepsize() const { return (span_[1] - span_[0]) / N_; }
};
namespace PhasePressODE {
template <class Density>
class Water {
public:
Water(const Density& rho,
const int np,
const int ix,
const double norm_grav)
: rho_(rho)
, svol_(np, 0)
, ix_(ix)
, g_(norm_grav)
{
svol_[ix_] = 1.0;
}
double
operator()(const double /* depth */,
const double press) const
{
return this->density(press) * g_;
}
private:
const Density& rho_;
std::vector<double> svol_;
const int ix_;
const double g_;
double
density(const double press) const
{
const std::vector<double>& rho = rho_(press, svol_);
return rho[ix_];
}
};
template <class Density, class RS>
class Oil {
public:
Oil(const Density& rho,
const RS& rs,
const int np,
const int oix,
const int gix,
const double norm_grav)
: rho_(rho)
, rs_(rs)
, svol_(np, 0)
, oix_(oix)
, gix_(gix)
, g_(norm_grav)
{
svol_[oix_] = 1.0;
}
double
operator()(const double depth,
const double press) const
{
return this->density(depth, press) * g_;
}
private:
const Density& rho_;
const RS& rs_;
mutable std::vector<double> svol_;
const int oix_;
const int gix_;
const double g_;
double
density(const double depth,
const double press) const
{
if (gix_ >= 0) {
svol_[gix_] = rs_(depth, press);
}
const std::vector<double>& rho = rho_(press, svol_);
return rho[oix_];
}
};
template <class Density, class RV>
class Gas {
public:
Gas(const Density& rho,
const RV& rv,
const int np,
const int gix,
const int oix,
const double norm_grav)
: rho_(rho)
, rv_(rv)
, svol_(np, 0)
, gix_(gix)
, oix_(oix)
, g_(norm_grav)
{
svol_[gix_] = 1.0;
}
double
operator()(const double depth,
const double press) const
{
return this->density(depth, press) * g_;
}
private:
const Density& rho_;
const RV& rv_;
mutable std::vector<double> svol_;
const int gix_;
const int oix_;
const double g_;
double
density(const double depth,
const double press) const
{
if (oix_ >= 0) {
svol_[oix_] = rv_(depth, press);
}
const std::vector<double>& rho = rho_(press, svol_);
return rho[gix_];
}
};
} // namespace PhasePressODE
namespace PhaseUsed {
inline bool
water(const PhaseUsage& pu)
{
return bool(pu.phase_used[ Opm::BlackoilPhases::Aqua ]);
}
inline bool
oil(const PhaseUsage& pu)
{
return bool(pu.phase_used[ Opm::BlackoilPhases::Liquid ]);
}
inline bool
gas(const PhaseUsage& pu)
{
return bool(pu.phase_used[ Opm::BlackoilPhases::Vapour ]);
}
} // namespace PhaseUsed
namespace PhaseIndex {
inline int
water(const PhaseUsage& pu)
{
int i = -1;
if (PhaseUsed::water(pu)) {
i = pu.phase_pos[ Opm::BlackoilPhases::Aqua ];
}
return i;
}
inline int
oil(const PhaseUsage& pu)
{
int i = -1;
if (PhaseUsed::oil(pu)) {
i = pu.phase_pos[ Opm::BlackoilPhases::Liquid ];
}
return i;
}
inline int
gas(const PhaseUsage& pu)
{
int i = -1;
if (PhaseUsed::gas(pu)) {
i = pu.phase_pos[ Opm::BlackoilPhases::Vapour ];
}
return i;
}
} // namespace PhaseIndex
namespace PhasePressure {
template <class PressFunction,
class CellRange>
void
assign(const UnstructuredGrid& G ,
const std::array<PressFunction, 2>& f ,
const double split,
const CellRange& cells,
std::vector<double>& p )
{
const int nd = G.dimensions;
enum { up = 0, down = 1 };
std::vector<double>::size_type c = 0;
for (typename CellRange::const_iterator
ci = cells.begin(), ce = cells.end();
ci != ce; ++ci, ++c)
{
assert (c < p.size());
const double z = G.cell_centroids[(*ci)*nd + (nd - 1)];
p[c] = (z < split) ? f[up](z) : f[down](z);
}
}
template <class Region,
class CellRange>
void
water(const UnstructuredGrid& G ,
const Region& reg ,
const std::array<double,2>& span ,
const double grav ,
const double po_woc,
const CellRange& cells ,
std::vector<double>& press )
{
using PhasePressODE::Water;
typedef Water<typename Region::CalcDensity> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int wix = PhaseIndex::water(pu);
ODE drho(reg.densityCalculator(), pu.num_phases, wix, grav);
const double z0 = reg.zwoc();
const double p0 = po_woc - reg.pcow_woc(); // Pcow = Po - Pw
std::array<double,2> up = {{ z0, span[0] }};
std::array<double,2> down = {{ z0, span[1] }};
typedef Details::RK4IVP<ODE> WPress;
std::array<WPress,2> wpress = {
{
WPress(drho, up , p0, 100)
,
WPress(drho, down, p0, 100)
}
};
assign(G, wpress, z0, cells, press);
}
template <class Region,
class CellRange>
void
oil(const UnstructuredGrid& G ,
const Region& reg ,
const std::array<double,2>& span ,
const double grav ,
const CellRange& cells ,
std::vector<double>& press ,
double& po_woc,
double& po_goc)
{
using PhasePressODE::Oil;
typedef Oil<typename Region::CalcDensity,
typename Region::CalcDissolution> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int oix = PhaseIndex::oil(pu);
const int gix = PhaseIndex::gas(pu);
ODE drho(reg.densityCalculator(),
reg.dissolutionCalculator(),
pu.num_phases, oix, gix, grav);
const double z0 = reg.datum();
const double p0 = reg.pressure();
std::array<double,2> up = {{ z0, span[0] }};
std::array<double,2> down = {{ z0, span[1] }};
typedef Details::RK4IVP<ODE> OPress;
std::array<OPress,2> opress = {
{
OPress(drho, up , p0, 100)
,
OPress(drho, down, p0, 100)
}
};
assign(G, opress, z0, cells, press);
po_woc = -std::numeric_limits<double>::max();
po_goc = -std::numeric_limits<double>::max();
const double woc = reg.zwoc();
// Compute Oil pressure at WOC
if (z0 > woc) { po_woc = opress[0](woc); } // WOC above datum
else if (z0 < woc) { po_woc = opress[1](woc); } // WOC below datum
else { po_woc = p0; } // WOC *at* datum
const double goc = reg.zgoc();
// Compute Oil pressure at GOC
if (z0 > goc) { po_goc = opress[0](goc); } // GOC above datum
else if (z0 < goc) { po_goc = opress[1](goc); } // GOC below datum
else { po_goc = p0; } // GOC *at* datum
}
template <class Region,
class CellRange>
void
gas(const UnstructuredGrid& G ,
const Region& reg ,
const std::array<double,2>& span ,
const double grav ,
const double po_goc,
const CellRange& cells ,
std::vector<double>& press )
{
using PhasePressODE::Gas;
typedef Gas<typename Region::CalcDensity,
typename Region::CalcEvaporation> ODE;
const PhaseUsage& pu = reg.phaseUsage();
const int gix = PhaseIndex::gas(pu);
const int oix = PhaseIndex::oil(pu);
ODE drho(reg.densityCalculator(),
reg.evaporationCalculator(),
pu.num_phases, gix, oix, grav);
const double z0 = reg.zgoc();
const double p0 = po_goc + reg.pcgo_goc(); // Pcog = Pg - Po
std::array<double,2> up = {{ z0, span[0] }};
std::array<double,2> down = {{ z0, span[1] }};
typedef Details::RK4IVP<ODE> GPress;
std::array<GPress,2> gpress = {
{
GPress(drho, up , p0, 100)
,
GPress(drho, down, p0, 100)
}
};
assign(G, gpress, z0, cells, press);
}
} // namespace PhasePressure
template <class Region,
class CellRange>
void
equilibrateOWG(const UnstructuredGrid& G,
const Region& reg,
const double grav,
const std::array<double,2>& span,
const CellRange& cells,
std::vector< std::vector<double> >& press)
{
const PhaseUsage& pu = reg.phaseUsage();
double po_woc = -1, po_goc = -1;
if (PhaseUsed::oil(pu)) {
const int oix = PhaseIndex::oil(pu);
PhasePressure::oil(G, reg, span, grav, cells,
press[ oix ], po_woc, po_goc);
}
if (PhaseUsed::water(pu)) {
const int wix = PhaseIndex::water(pu);
PhasePressure::water(G, reg, span, grav, po_woc,
cells, press[ wix ]);
}
if (PhaseUsed::gas(pu)) {
const int gix = PhaseIndex::gas(pu);
PhasePressure::gas(G, reg, span, grav, po_goc,
cells, press[ gix ]);
}
}
} // namespace Details
namespace Equil {
template <class Region,
class CellRange>
std::vector< std::vector<double> >
phasePressures(const UnstructuredGrid& G,
const Region& reg,
const CellRange& cells,
const double grav)
{
std::array<double,2> span =
{{ std::numeric_limits<double>::max() ,
-std::numeric_limits<double>::max() }}; // Symm. about 0.
int ncell = 0;
{
// This code is only supported in three space dimensions
assert (G.dimensions == 3);
const int nd = G.dimensions;
// Define short-name aliases to reduce visual clutter.
const double* const nc = & G.node_coordinates[0];
const int* const cfp = & G.cell_facepos[0];
const int* const cf = & G.cell_faces[0];
const int* const fnp = & G.face_nodepos[0];
const int* const fn = & G.face_nodes[0];
// Define vertical span as
//
// [minimum(node depth(cells)), maximum(node depth(cells))]
//
// Note: We use a sledgehammer approach--looping all
// the nodes of all the faces of all the 'cells'--to
// compute those bounds. This necessarily entails
// visiting some nodes (and faces) multiple times.
//
// Note: The implementation of 'RK4IVP<>' implicitly
// imposes the requirement that cell centroids are all
// within this vertical span. That requirement is not
// checked.
for (typename CellRange::const_iterator
ci = cells.begin(), ce = cells.end();
ci != ce; ++ci, ++ncell)
{
for (const int
*fi = & cf[ cfp[*ci + 0] ],
*fe = & cf[ cfp[*ci + 1] ];
fi != fe; ++fi)
{
for (const int
*i = & fn[ fnp[*fi + 0] ],
*e = & fn[ fnp[*fi + 1] ];
i != e; ++i)
{
const double z = nc[(*i)*nd + (nd - 1)];
if (z < span[0]) { span[0] = z; }
if (z > span[1]) { span[1] = z; }
}
}
}
}
const int np = reg.phaseUsage().num_phases;
typedef std::vector<double> pval;
std::vector<pval> press(np, pval(ncell, 0.0));
const double z0 = reg.datum();
const double zwoc = reg.zwoc ();
const double zgoc = reg.zgoc ();
if (! ((zgoc > z0) || (z0 > zwoc))) {
// Datum depth in oil zone (zgoc <= z0 <= zwoc)
Details::equilibrateOWG(G, reg, grav, span, cells, press);
} else {
OPM_THROW(std::runtime_error, "Cannot initialise: the datum depth must be in the oil zone.");
}
return press;
}
template <class Region, class CellRange>
std::vector< std::vector<double> >
phaseSaturations(const Region& reg,
const CellRange& cells,
const BlackoilPropertiesInterface& props,
std::vector< std::vector<double> >& phase_pressures)
{
const double z0 = reg.datum();
const double zwoc = reg.zwoc ();
const double zgoc = reg.zgoc ();
if ((zgoc > z0) || (z0 > zwoc)) {
OPM_THROW(std::runtime_error, "Cannot initialise: the datum depth must be in the oil zone.");
}
if (!reg.phaseUsage().phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Cannot initialise: not handling water-gas cases.");
}
std::vector< std::vector<double> > phase_saturations = phase_pressures; // Just to get the right size.
double smin[BlackoilPhases::MaxNumPhases] = { 0.0 };
double smax[BlackoilPhases::MaxNumPhases] = { 0.0 };
const bool water = reg.phaseUsage().phase_used[BlackoilPhases::Aqua];
const bool gas = reg.phaseUsage().phase_used[BlackoilPhases::Vapour];
const int oilpos = reg.phaseUsage().phase_pos[BlackoilPhases::Liquid];
const int waterpos = reg.phaseUsage().phase_pos[BlackoilPhases::Aqua];
const int gaspos = reg.phaseUsage().phase_pos[BlackoilPhases::Vapour];
std::vector<double>::size_type local_index = 0;
for (typename CellRange::const_iterator ci = cells.begin(); ci != cells.end(); ++ci, ++local_index) {
const int cell = *ci;
props.satRange(1, &cell, smin, smax);
// Find saturations from pressure differences by
// inverting capillary pressure functions.
double sw = 0.0;
if (water) {
const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
sw = satFromPc(props, waterpos, cell, pcov);
phase_saturations[waterpos][local_index] = sw;
}
double sg = 0.0;
if (gas) {
// Note that pcog is defined to be (pg - po), not (po - pg).
const double pcog = phase_pressures[gaspos][local_index] - phase_pressures[oilpos][local_index];
const double increasing = true; // pcog(sg) expected to be increasing function
sg = satFromPc(props, gaspos, cell, pcog, increasing);
phase_saturations[gaspos][local_index] = sg;
}
bool overlap = false;
if (gas && water && (sg + sw > 1.0)) {
// Overlapping gas-oil and oil-water transition
// zones can lead to unphysical saturations when
// treated as above. Must recalculate using gas-water
// capillary pressure.
const double pcgw = phase_pressures[gaspos][local_index] - phase_pressures[waterpos][local_index];
sw = satFromSumOfPcs(props, waterpos, gaspos, cell, pcgw);
sg = 1.0 - sw;
phase_saturations[waterpos][local_index] = sw;
phase_saturations[gaspos][local_index] = sg;
overlap = true;
}
phase_saturations[oilpos][local_index] = 1.0 - sw - sg;
// Adjust phase pressures for max and min saturation ...
double pc[BlackoilPhases::MaxNumPhases];
double sat[BlackoilPhases::MaxNumPhases];
if (sw > smax[waterpos]-1.0e-6) {
sat[waterpos] = smax[waterpos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[oilpos][local_index] = phase_pressures[waterpos][local_index] + pc[waterpos];
} else if (overlap || sg > smax[gaspos]-1.0e-6) {
sat[gaspos] = smax[gaspos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[oilpos][local_index] = phase_pressures[gaspos][local_index] - pc[gaspos];
}
if (sg < smin[gaspos]+1.0e-6) {
sat[gaspos] = smin[gaspos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[gaspos][local_index] = phase_pressures[oilpos][local_index] + pc[gaspos];
}
if (sw < smin[waterpos]+1.0e-6) {
sat[waterpos] = smin[waterpos];
props.capPress(1, sat, &cell, pc, 0);
phase_pressures[waterpos][local_index] = phase_pressures[oilpos][local_index] - pc[waterpos];
}
}
return phase_saturations;
}
/**
* Compute initial Rs values.
*
* \tparam CellRangeType Type of cell range that demarcates the
* cells pertaining to the current
* equilibration region. Must implement
* methods begin() and end() to bound the range
* as well as provide an inner type,
* const_iterator, to traverse the range.
*
* \param[in] grid Grid.
* \param[in] cells Range that spans the cells of the current
* equilibration region.
* \param[in] oil_pressure Oil pressure for each cell in range.
* \param[in] rs_func Rs as function of pressure and depth.
* \return Rs values, one for each cell in the 'cells' range.
*/
template <class CellRangeType>
std::vector<double> computeRs(const UnstructuredGrid& grid,
const CellRangeType& cells,
const std::vector<double> oil_pressure,
const Miscibility::RsFunction& rs_func,
const std::vector<double> gas_saturation)
{
assert(grid.dimensions == 3);
std::vector<double> rs(cells.size());
int count = 0;
for (auto it = cells.begin(); it != cells.end(); ++it, ++count) {
const double depth = grid.cell_centroids[3*(*it) + 2];
rs[count] = rs_func(depth, oil_pressure[count], gas_saturation[count]);
}
return rs;
}
} // namespace Equil
namespace
{
/// Convert saturations from a vector of individual phase saturation vectors
/// to an interleaved format where all values for a given cell come before all
/// values for the next cell, all in a single vector.
std::vector<double> convertSats(const std::vector< std::vector<double> >& sat)
{
const int np = sat.size();
const int nc = sat[0].size();
std::vector<double> s(np * nc);
for (int c = 0; c < nc; ++c) {
for (int p = 0; p < np; ++p) {
s[np*c + p] = sat[p][c];
}
}
return s;
}
}
/**
* Compute initial state by an equilibration procedure.
*
* The following state fields are modified:
* pressure(),
* saturation(),
* surfacevol(),
* gasoilratio(),
* rv().
*
* \param[in] grid Grid.
* \param[in] props Property object, pvt and capillary properties are used.
* \param[in] deck Simulation deck, used to obtain EQUIL and related data.
* \param[in] gravity Acceleration of gravity, assumed to be in Z direction.
*/
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const EclipseGridParser& deck,
const double gravity,
BlackoilState& state)
{
typedef Equil::DeckDependent::InitialStateComputer<EclipseGridParser> ISC;
ISC isc(props, deck, grid, gravity);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua];
state.pressure() = isc.press()[ref_phase];
state.saturation() = convertSats(isc.saturation());
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// TODO: state.surfacevol() must be computed from s, rs, rv.
}
void initStateEquil(const UnstructuredGrid& grid,
const BlackoilPropertiesInterface& props,
const Opm::DeckConstPtr newParserDeck,
const double gravity,
BlackoilState& state)
{
typedef Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> ISC;
ISC isc(props, newParserDeck, grid, gravity);
const auto pu = props.phaseUsage();
const int ref_phase = pu.phase_used[BlackoilPhases::Liquid]
? pu.phase_pos[BlackoilPhases::Liquid]
: pu.phase_pos[BlackoilPhases::Aqua];
state.pressure() = isc.press()[ref_phase];
state.saturation() = convertSats(isc.saturation());
state.gasoilratio() = isc.rs();
state.rv() = isc.rv();
// TODO: state.surfacevol() must be computed from s, rs, rv.
}
} // namespace Opm
#endif // OPM_INITSTATEEQUIL_IMPL_HEADER_INCLUDED

205
opm/core/utility/RegionMapping.hpp Normal file
View File

@ -0,0 +1,205 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef OPM_REGIONMAPPING_HEADER_INCLUDED
#define OPM_REGIONMAPPING_HEADER_INCLUDED
#include <vector>
namespace Opm
{
/**
* Forward and reverse mappings between cells and
* regions/partitions (e.g., the ECLIPSE-style 'SATNUM',
* 'PVTNUM', or 'EQUILNUM' arrays).
*
* \tparam Region Type of a forward region mapping. Expected
* to provide indexed access through
* operator[]() as well as inner types
* 'value_type', 'size_type', and
* 'const_iterator'.
*/
template < class Region = std::vector<int> >
class RegionMapping {
public:
/**
* Constructor.
*
* \param[in] reg Forward region mapping, restricted to
* active cells only.
*/
explicit
RegionMapping(const Region& reg)
: reg_(reg)
{
rev_.init(reg_);
}
/**
* Type of forward (cell-to-region) mapping result.
* Expected to be an integer.
*/
typedef typename Region::value_type RegionId;
/**
* Type of reverse (region-to-cell) mapping (element)
* result.
*/
typedef typename Region::size_type CellId;
/**
* Type of reverse region-to-cell range bounds and
* iterators.
*/
typedef typename std::vector<CellId>::const_iterator CellIter;
/**
* Range of cells. Result from reverse (region-to-cell)
* mapping.
*/
class CellRange {
public:
/**
* Constructor.
*
* \param[in] b Beginning of range.
* \param[in] e One past end of range.
*/
CellRange(const CellIter b,
const CellIter e)
: b_(b), e_(e)
{}
/**
* Read-only iterator on cell ranges.
*/
typedef CellIter const_iterator;
/**
* Size type for this range.
*/
typedef typename std::vector<CellId>::size_type size_type;
/**
* Beginning of cell range.
*/
const_iterator begin() const { return b_; }
/**
* One past end of cell range.
*/
const_iterator end() const { return e_; }
/**
* Number of elements in the range.
*/
size_type size() const { return e_ - b_; }
private:
const_iterator b_;
const_iterator e_;
};
/**
* Number of declared regions in cell-to-region mapping.
*/
RegionId
numRegions() const { return RegionId(rev_.p.size()) - 1; }
/**
* Compute region number of given active cell.
*
* \param[in] c Active cell
* \return Region to which @c c belongs.
*/
RegionId
region(const CellId c) const { return reg_[c]; }
/**
* Extract active cells in particular region.
*
* \param[in] r Region number
* \returns Range of active cells in region @c r.
*/
CellRange
cells(const RegionId r) const {
const RegionId i = r - rev_.low;
return CellRange(rev_.c.begin() + rev_.p[i + 0],
rev_.c.begin() + rev_.p[i + 1]);
}
private:
/**
* Copy of forward region mapping (cell-to-region).
*/
Region reg_;
/**
* Reverse mapping (region-to-cell).
*/
struct {
typedef typename std::vector<CellId>::size_type Pos;
std::vector<Pos> p; /**< Region start pointers */
std::vector<CellId> c; /**< Region cells */
RegionId low; /**< Smallest region number */
/**
* Compute reverse mapping. Standard linear insertion
* sort algorithm.
*/
void
init(const Region& reg)
{
typedef typename Region::const_iterator CI;
const std::pair<CI,CI>
m = std::minmax_element(reg.begin(), reg.end());
low = *m.first;
const typename Region::size_type
n = *m.second - low + 1;
p.resize(n + 1); std::fill(p.begin(), p.end(), Pos(0));
for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
p[ reg[i] - low + 1 ] += 1;
}
for (typename std::vector<Pos>::size_type
i = 1, sz = p.size(); i < sz; ++i) {
p[0] += p[i];
p[i] = p[0] - p[i];
}
assert (p[0] ==
static_cast<typename Region::size_type>(reg.size()));
c.resize(reg.size());
for (CellId i = 0, nc = reg.size(); i < nc; ++i) {
c[ p[ reg[i] - low + 1 ] ++ ] = i;
}
p[0] = 0;
}
} rev_; /**< Reverse mapping instance */
};
} // namespace Opm
#endif // OPM_REGIONMAPPING_HEADER_INCLUDED

6
opm/core/wells.h
View File

@ -65,8 +65,10 @@ struct Wells
/**
* Component fractions for each well. Array of size
* <CODE>number_of_wells * number_of_phases</CODE>.
* This is intended to be used for injection wells. For production wells
* the component fractions will vary and cannot be specified a priori.
* For injection wells, this gives the injected component mix.
* For production wells the component fractions of the wellbore
* will vary and cannot be specified a priori, the component mix
* given here should be considered a default or preferred mix.
*/
double *comp_frac;

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

@ -550,6 +550,30 @@ namespace Opm
OPM_THROW(std::runtime_error, "Control mode type " << mode << " not present in well " << well_names[well_index]);
}
set_current_control(well_index, cpos, w_);
// Set well component fraction to match preferred phase for the well.
double cf[3] = { 0.0, 0.0, 0.0 };
{
switch (well->getPreferredPhase()) {
case Phase::WATER:
if (!phaseUsage.phase_used[BlackoilPhases::Aqua]) {
OPM_THROW(std::runtime_error, "Water phase not used, yet found water-preferring well.");
}
cf[phaseUsage.phase_pos[BlackoilPhases::Aqua]] = 1.0;
break;
case Phase::OIL:
if (!phaseUsage.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Oil phase not used, yet found oil-preferring well.");
}
cf[phaseUsage.phase_pos[BlackoilPhases::Liquid]] = 1.0;
case Phase::GAS:
if (!phaseUsage.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "Gas phase not used, yet found gas-preferring well.");
}
cf[phaseUsage.phase_pos[BlackoilPhases::Vapour]] = 1.0;
}
std::copy(cf, cf + phaseUsage.num_phases, w_->comp_frac + well_index*phaseUsage.num_phases);
}
}
well_index++;
}

38
tests/capillary.DATA Normal file
View File

@ -0,0 +1,38 @@
WATER
OIL
GAS
TABDIMS
1 1 40 20 1 20 /
EQLDIMS
-- NTEQUL
1 /
PVDO
100 1.0 1.0
200 0.9 1.0
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
SWOF
0.2 0 1 0.4
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
DENSITY
700 1000 1
/
EQUIL
50 150 50 0.25 20 0.35 1* 1* 0
/

128
tests/capillary_overlap.DATA Normal file
View File

@ -0,0 +1,128 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
PORO
20*0.2
/
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVDO
100 1.0 1.0
200 0.9 1.0
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1* 1* 0
/
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RESTART=3' 'NEWTON=2' /
END

38
tests/deadfluids.DATA Normal file
View File

@ -0,0 +1,38 @@
WATER
OIL
GAS
TABDIMS
1 1 40 20 1 20 /
EQLDIMS
-- NTEQUL
1 /
PVDO
100 1.0 1.0
200 0.5 1.0
/
PVDG
100 0.05 0.1
200 0.02 0.2
/
SWOF
0 0 1 0
1 1 0 0
/
SGOF
0 0 1 0
1 1 0 0
/
DENSITY
700 1000 10
/
EQUIL
5 150 5 0 2 0 1* 1* 0
/

131
tests/equil_livegas.DATA Normal file
View File

@ -0,0 +1,131 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
PORO
20*0.2
/
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVDO
100 1.0 1.0
200 0.9 1.0
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1* 1* 0
/
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

140
tests/equil_liveoil.DATA Normal file
View File

@ -0,0 +1,140 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
PORO
20*0.2
/
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1* 1* 0
/
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

65
tests/equil_liveoil_grid.DATA Normal file
View File

@ -0,0 +1,65 @@
WATER
OIL
GAS
DISGAS
DIMENS
1 1 20
/
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
TOPS
4*0.0
/
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
100 40. 1.0120 1.17 /
200 80. 1.0255 1.14 /
300 120. 1.0380 1.11 /
400 160. 1.0510 1.08 /
500 200. 1.0630 1.06 /
600 240. 1.0750 1.03 /
700 280. 1.0870 1.00 /
800 320. 1.0985 .98 /
900 360. 1.1100 .95 /
1000 400. 1.1200 .94
500. 1.1189 .94 /
/
PVDG
100 0.010 0.1
200 0.005 0.2
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
DENSITY
700 1000 1
/
EQUIL
45 150 50 0.25 45 0.35
/

151
tests/equil_rsvd_and_rvvd.DATA Normal file
View File

@ -0,0 +1,151 @@
NOECHO
RUNSPEC ======
WATER
OIL
GAS
DISGAS
VAPOIL
TABDIMS
1 1 40 20 1 20 /
DIMENS
1 1 20
/
WELLDIMS
30 10 2 30 /
START
1 'JAN' 1990 /
NSTACK
25 /
EQLDIMS
-- NTEQUL
1 /
FMTOUT
FMTIN
GRID ======
DXV
1.0
/
DYV
1.0
/
DZV
20*5.0
/
PORO
20*0.2
/
PERMZ
20*1.0
/
PERMY
20*100.0
/
PERMX
20*100.0
/
BOX
1 1 1 1 1 1 /
TOPS
0.0
/
PROPS ======
PVTO
-- Rs Pbub Bo Vo
0 1. 1.0000 1.20 /
20 40. 1.0120 1.17 /
40 80. 1.0255 1.14 /
60 120. 1.0380 1.11 /
80 160. 1.0510 1.08 /
100 200. 1.0630 1.06 /
120 240. 1.0750 1.03 /
140 280. 1.0870 1.00 /
160 320. 1.0985 .98 /
180 360. 1.1100 .95 /
200 400. 1.1200 .94
500. 1.1189 .94 /
/
PVTG
-- Pg Rv Bg Vg
100 0.0001 0.010 0.1
0.0 0.0104 0.1 /
200 0.0004 0.005 0.2
0.0 0.0054 0.2 /
/
SWOF
0.2 0 1 0.9
1 1 0 0.1
/
SGOF
0 0 1 0.2
0.8 1 0 0.5
/
PVTW
--RefPres Bw Comp Vw Cv
1. 1.0 4.0E-5 0.96 0.0 /
ROCK
--RefPres Comp
1. 5.0E-5 /
DENSITY
700 1000 1
/
SOLUTION ======
EQUIL
45 150 50 0.25 45 0.35 1 1 0
/
RSVD
0 0.0
100 100. /
RVVD
0. 0.
100. 0.0001 /
RPTSOL
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=2' /
SUMMARY ======
RUNSUM
SEPARATE
SCHEDULE ======
RPTSCHED
'PRES' 'PGAS' 'PWAT' 'SOIL' 'SWAT' 'SGAS' 'RS' 'RESTART=3' 'NEWTON=2' /
END

764
tests/test_equil.cpp Normal file
View File

@ -0,0 +1,764 @@
/*
Copyright 2014 SINTEF ICT, Applied Mathematics.
*/
#include "config.h"
/* --- Boost.Test boilerplate --- */
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE UnitsTest
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
/* --- our own headers --- */
#include <opm/core/simulator/initStateEquil.hpp>
#include <opm/core/grid.h>
#include <opm/core/grid/cart_grid.h>
#include <opm/core/grid/GridManager.hpp>
#include <opm/core/props/BlackoilPropertiesBasic.hpp>
#include <opm/core/props/BlackoilPropertiesFromDeck.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/core/pressure/msmfem/partition.h>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/core/utility/Units.hpp>
#include <array>
#include <iostream>
#include <limits>
#include <memory>
#include <numeric>
#include <sstream>
#include <string>
#include <vector>
BOOST_AUTO_TEST_SUITE ()
BOOST_AUTO_TEST_CASE (PhasePressure)
{
typedef std::vector<double> PVal;
typedef std::vector<PVal> PPress;
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::parameter::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::Equil::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
Opm::Equil::EquilRecord record =
{
{ 0 , 1e5 } , // Datum depth, pressure
{ 5 , 0 } , // Zwoc , Pcow_woc
{ 0 , 0 } // Zgoc , Pcgo_goc
};
Opm::Equil::EquilReg<RhoCalc>
region(record, calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage());
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);
const double grav = 10;
const PPress ppress = Opm::Equil::phasePressures(*G, region, cells, grav);
const int first = 0, last = G->number_of_cells - 1;
const double reltol = 1.0e-8;
BOOST_CHECK_CLOSE(ppress[0][first] , 90e3 , reltol);
BOOST_CHECK_CLOSE(ppress[0][last ] , 180e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
}
BOOST_AUTO_TEST_CASE (CellSubset)
{
typedef std::vector<double> PVal;
typedef std::vector<PVal> PPress;
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::parameter::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::Equil::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
Opm::Equil::EquilRecord record[] =
{
{
{ 0 , 1e5 } , // Datum depth, pressure
{ 2.5 , -0.075e5 } , // Zwoc , Pcow_woc
{ 0 , 0 } // Zgoc , Pcgo_goc
}
,
{
{ 5 , 1.35e5 } , // Datum depth, pressure
{ 7.5 , -0.225e5 } , // Zwoc , Pcow_woc
{ 5 , 0 } // Zgoc , Pcgo_goc
}
};
Opm::Equil::EquilReg<RhoCalc> region[] =
{
Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
};
const int cdim[] = { 2, 1, 2 };
int ncoarse = cdim[0];
for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
std::vector< std::vector<int> > cells(ncoarse);
for (int c = 0; c < G->number_of_cells; ++c) {
int ci = c;
const int i = ci % G->cartdims[0]; ci /= G->cartdims[0];
const int j = ci % G->cartdims[1];
const int k = ci / G->cartdims[1];
const int ic = (i / (G->cartdims[0] / cdim[0]));
const int jc = (j / (G->cartdims[1] / cdim[1]));
const int kc = (k / (G->cartdims[2] / cdim[2]));
const int ix = ic + cdim[0]*(jc + cdim[1]*kc);
assert ((0 <= ix) && (ix < ncoarse));
cells[ix].push_back(c);
}
PPress ppress(2, PVal(G->number_of_cells, 0));
for (std::vector< std::vector<int> >::const_iterator
r = cells.begin(), e = cells.end();
r != e; ++r)
{
const int rno = int(r - cells.begin());
const double grav = 10;
const PPress p =
Opm::Equil::phasePressures(*G, region[rno], *r, grav);
PVal::size_type i = 0;
for (std::vector<int>::const_iterator
c = r->begin(), ce = r->end();
c != ce; ++c, ++i)
{
assert (i < p[0].size());
ppress[0][*c] = p[0][i];
ppress[1][*c] = p[1][i];
}
}
const int first = 0, last = G->number_of_cells - 1;
const double reltol = 1.0e-8;
BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol);
BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
}
BOOST_AUTO_TEST_CASE (RegMapping)
{
typedef std::vector<double> PVal;
typedef std::vector<PVal> PPress;
std::shared_ptr<UnstructuredGrid>
G(create_grid_cart3d(10, 1, 10), destroy_grid);
Opm::parameter::ParameterGroup param;
{
using Opm::unit::kilogram;
using Opm::unit::meter;
using Opm::unit::cubic;
std::stringstream dens; dens << 700*kilogram/cubic(meter);
param.insertParameter("rho2", dens.str());
}
typedef Opm::BlackoilPropertiesBasic Props;
Props props(param, G->dimensions, G->number_of_cells);
typedef Opm::Equil::DensityCalculator<Opm::BlackoilPropertiesInterface> RhoCalc;
RhoCalc calc(props, 0);
Opm::Equil::EquilRecord record[] =
{
{
{ 0 , 1e5 } , // Datum depth, pressure
{ 2.5 , -0.075e5 } , // Zwoc , Pcow_woc
{ 0 , 0 } // Zgoc , Pcgo_goc
}
,
{
{ 5 , 1.35e5 } , // Datum depth, pressure
{ 7.5 , -0.225e5 } , // Zwoc , Pcow_woc
{ 5 , 0 } // Zgoc , Pcgo_goc
}
};
Opm::Equil::EquilReg<RhoCalc> region[] =
{
Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::Equil::EquilReg<RhoCalc>(record[0], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
,
Opm::Equil::EquilReg<RhoCalc>(record[1], calc,
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
std::make_shared<Opm::Equil::Miscibility::NoMixing>(),
props.phaseUsage())
};
std::vector<int> eqlnum(G->number_of_cells);
{
std::vector<int> cells(G->number_of_cells);
std::iota(cells.begin(), cells.end(), 0);
const int cdim[] = { 2, 1, 2 };
int ncoarse = cdim[0];
for (std::size_t d = 1; d < 3; ++d) { ncoarse *= cdim[d]; }
partition_unif_idx(G->dimensions, G->number_of_cells,
G->cartdims, cdim,
&cells[0], &eqlnum[0]);
}
Opm::RegionMapping<> eqlmap(eqlnum);
PPress ppress(2, PVal(G->number_of_cells, 0));
for (int r = 0, e = eqlmap.numRegions(); r != e; ++r)
{
const Opm::RegionMapping<>::CellRange&
rng = eqlmap.cells(r);
const int rno = r;
const double grav = 10;
const PPress p =
Opm::Equil::phasePressures(*G, region[rno], rng, grav);
PVal::size_type i = 0;
for (Opm::RegionMapping<>::CellRange::const_iterator
c = rng.begin(), ce = rng.end();
c != ce; ++c, ++i)
{
assert (i < p[0].size());
ppress[0][*c] = p[0][i];
ppress[1][*c] = p[1][i];
}
}
const int first = 0, last = G->number_of_cells - 1;
const double reltol = 1.0e-8;
BOOST_CHECK_CLOSE(ppress[0][first] , 105e3 , reltol);
BOOST_CHECK_CLOSE(ppress[0][last ] , 195e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][first] , 103.5e3 , reltol);
BOOST_CHECK_CLOSE(ppress[1][last ] , 166.5e3 , reltol);
}
BOOST_AUTO_TEST_CASE (DeckAllDead)
{
std::shared_ptr<UnstructuredGrid>
grid(create_grid_cart3d(1, 1, 10), destroy_grid);
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("deadfluids.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, *grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, *grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid->number_of_cells);
const int first = 0, last = grid->number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-3;
BOOST_CHECK_CLOSE(pressures[0][first] , 1.496329839e7 , reltol);
BOOST_CHECK_CLOSE(pressures[0][last ] , 1.50473245e7 , reltol);
BOOST_CHECK_CLOSE(pressures[1][last] , 1.50473245e7 , reltol);
}
BOOST_AUTO_TEST_CASE (CapillaryInversion)
{
// Test setup.
Opm::GridManager gm(1, 1, 40, 1.0, 1.0, 2.5);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
// Test the capillary inversion for oil-water.
const int cell = 0;
const double reltol = 1.0e-7;
{
const int phase = 0;
const bool increasing = false;
const std::vector<double> pc = { 10.0e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.099e5, 0.0e5, -10.0e5 };
const std::vector<double> s = { 0.2, 0.2, 0.2, 0.466666666666, 0.733333333333, 1.0, 1.0, 1.0, 1.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::Equil::satFromPc(props, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
// Test the capillary inversion for gas-oil.
{
const int phase = 2;
const bool increasing = true;
const std::vector<double> pc = { 10.0e5, 0.6e5, 0.5e5, 0.4e5, 0.3e5, 0.2e5, 0.1e5, 0.0e5, -10.0e5 };
const std::vector<double> s = { 0.8, 0.8, 0.8, 0.533333333333, 0.266666666666, 0.0, 0.0, 0.0, 0.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::Equil::satFromPc(props, phase, cell, pc[i], increasing);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
// Test the capillary inversion for gas-water.
{
const int water = 0;
const int gas = 2;
const std::vector<double> pc = { 0.9e5, 0.8e5, 0.6e5, 0.4e5, 0.3e5 };
const std::vector<double> s = { 0.2, 0.333333333333, 0.6, 0.866666666666, 1.0 };
BOOST_REQUIRE(pc.size() == s.size());
for (size_t i = 0; i < pc.size(); ++i) {
const double s_computed = Opm::Equil::satFromSumOfPcs(props, water, gas, cell, pc[i]);
BOOST_CHECK_CLOSE(s_computed, s[i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (DeckWithCapillary)
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("capillary.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 10.0);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
BOOST_CHECK_CLOSE(pressures[0][first] , 1.469769063e7 , reltol);
BOOST_CHECK_CLOSE(pressures[0][last ] , 1.545e7 , reltol);
BOOST_CHECK_CLOSE(pressures[1][last] , 1.546e7 , reltol);
const auto& sats = comp.saturation();
const std::vector<double> s[3]{
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.425893333333, 0.774026666666, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0.00736, 0.792746666666, 0.8, 0.8, 0.8, 0.8, 0.574106666666, 0.225973333333, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.79264, 0.007253333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE(sats[phase].size() == s[phase].size());
for (size_t i = 0; i < s[phase].size(); ++i) {
BOOST_CHECK_CLOSE(sats[phase][i], s[phase][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (DeckWithCapillaryOverlap)
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("capillary_overlap.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[0][first] , 14832467.14, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last ] , 15479883.47, reltol);
BOOST_CHECK_CLOSE(pressures[1][last ] , 15489883.47, reltol);
const auto& sats = comp.saturation();
// std::cout << "Saturations:\n";
// for (const auto& sat : sats) {
// for (const double s : sat) {
// std::cout << s << ' ';
// }
// std::cout << std::endl;
// }
const std::vector<double> s_ecl[3]{// eclipse
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22874042, 0.53397995, 0.78454906, 0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20039, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77125955, 0.46602005, 0.015063271, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
const std::vector<double> s_opm[3]{ // opm
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2289309090909091, 0.53406545454545451, 0.78458, 0.9154, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2002466666666666, 0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77106909090909093, 0.46593454545454549, 0.015173333333333336, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
}
}
}
BOOST_AUTO_TEST_CASE (DeckWithLiveOil)
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("equil_liveoil.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48324e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49224e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[0][first], 1.483246714e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last], 1.547991652e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][first], 1.492246714e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][last], 1.548991652e7, reltol);
const auto& sats = comp.saturation();
// std::cout << "Saturations:\n";
// for (const auto& sat : sats) {
// for (const double s : sat) {
// std::cout << s << ' ';
// }
// std::cout << std::endl;
// }
const std::vector<double> s_ecl[3]{ // eclipse
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22898, 0.53422, 0.78470, 0.91531, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.20073, 0.08469, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77102, 0.46578, 0.01458, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
const std::vector<double> s_opm[3]{ // opm
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2291709091, 0.5343054545, 0.78472, 0.91529, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.2005866667, 0.08471, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.7708290909, 0.4656945455, 0.01469333333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], reltol);
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
}
std::cout << std::endl;
}
const auto& rs = comp.rs();
const std::vector<double> rs_opm {74.612335679539058, 74.649052116644228, 74.685786561426298, 74.722539022717172, // opm
74.759309509353145, 74.796098030174733, 74.8329045940269, 74.869729209758916,
74.906571886224327, 75.090675116639048, 75.0, 75.0,
75.0, 75.0, 75.0, 75.0,
75.0, 75.0, 75.0, 75.0};
const std::vector<double> rs_ecl {74.612228, 74.648956, 74.685707, 74.722473, // eclipse
74.759254, 74.796051, 74.832870, 74.875145,
74.969231, 75.090706, 75.000000, 75.000000,
75.000000, 75.000000, 75.000000, 75.000000,
75.000000, 75.000000, 75.000000, 75.000000};
for (size_t i = 0; i < rs_opm.size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(rs[i], rs_opm[i], reltol);
BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE (DeckWithLiveGas)
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("equil_livegas.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48215e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54801e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49115e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54901e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[0][first], 1.482150311e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last], 1.547988347e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][first], 1.491150311e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][last], 1.548988347e7, reltol);
const auto& sats = comp.saturation();
// std::cout << "Saturations:\n";
// for (const auto& sat : sats) {
// for (const double s : sat) {
// std::cout << s << ' ';
// }
// std::cout << std::endl;
// }
const std::vector<double> s_ecl[3]{ // eclipse
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24285614, 0.53869015, 0.78454906, 0.91542006, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18311, 0.08458, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75714386, 0.46130988, 0.032345835, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
const std::vector<double> s_opm[3]{ // opm
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.24310545, 0.5388, 0.78458, 0.91540, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.18288667, 0.0846, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.75689455, 0.4612, 0.03253333, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol);
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
}
std::cout << std::endl;
}
const auto& rv = comp.rv();
const std::vector<double> rv_opm { // opm
2.4884509e-4, 2.4910378e-4, 2.4936267e-4, 2.4962174e-4,
2.4988100e-4, 2.5014044e-4, 2.5040008e-4, 2.5065990e-4,
2.5091992e-4, 2.5118012e-4, 2.5223082e-4, 2.5105e-4,
2.5105e-4, 2.5105e-4, 2.5105e-4, 2.5105e-4,
2.5105e-4, 2.5105e-4, 2.5105e-4, 2.5105e-4};
const std::vector<double> rv_ecl { // eclipse
0.24884584E-03, 0.24910446E-03, 0.24936325E-03, 0.24962222E-03,
0.24988138E-03, 0.25014076E-03, 0.25040031E-03, 0.25066003E-03,
0.25091995E-03, 0.25118008E-03, 0.25223137E-03, 0.25104999E-03,
0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03,
0.25104999E-03, 0.25104999E-03, 0.25104999E-03, 0.25104999E-03};
for (size_t i = 0; i < rv_opm.size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], 100.*reltol);
BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_CASE (DeckWithRSVDAndRVVD)
{
Opm::GridManager gm(1, 1, 20, 1.0, 1.0, 5.0);
const UnstructuredGrid& grid = *(gm.c_grid());
Opm::ParserPtr parser(new Opm::Parser() );
Opm::DeckConstPtr deck = parser->parseFile("equil_rsvd_and_rvvd.DATA");
Opm::BlackoilPropertiesFromDeck props(deck, grid, false);
Opm::Equil::DeckDependent::InitialStateComputer<Opm::DeckConstPtr> comp(props, deck, grid, 9.80665);
const auto& pressures = comp.press();
BOOST_REQUIRE(pressures.size() == 3);
BOOST_REQUIRE(int(pressures[0].size()) == grid.number_of_cells);
const int first = 0, last = grid.number_of_cells - 1;
// The relative tolerance is too loose to be very useful,
// but the answer we are checking is the result of an ODE
// solver, and it is unclear if we should check it against
// the true answer or something else.
const double reltol = 1.0e-6;
const double reltol_ecl = 1.0;
BOOST_CHECK_CLOSE(pressures[0][first], 1.48350e+07, reltol_ecl); // eclipse
BOOST_CHECK_CLOSE(pressures[0][last], 1.54794e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][first], 1.49250e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[1][last], 1.54894e+07, reltol_ecl);
BOOST_CHECK_CLOSE(pressures[0][first], 1.483499660e7, reltol); // opm
BOOST_CHECK_CLOSE(pressures[0][last], 1.547924516e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][first], 1.492499660e7, reltol);
BOOST_CHECK_CLOSE(pressures[1][last], 1.548924516e7, reltol);
const auto& sats = comp.saturation();
// std::cout << "Saturations:\n";
// for (const auto& sat : sats) {
// for (const double s : sat) {
// std::cout << s << ' ';
// }
// std::cout << std::endl;
// }
const std::vector<double> s_ecl[3]{ // eclipse
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22206347, 0.52871972, 0.78150368, 0.91819441, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19656529, 0.081805572, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77793652, 0.47128031, 0.021931054, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
const std::vector<double> s_opm[3]{ // opm
{ 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.2, 0.22232000, 0.52882909, 0.78153000, 0.91817000, 1, 1, 1, 1, 1, 1, 1, 1, 1 },
{ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.19636333, 0.08183000, 0, 0, 0, 0, 0, 0, 0, 0, 0 },
{ 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.8, 0.77768000, 0.47117091, 0.02210667, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }
};
for (int phase = 0; phase < 3; ++phase) {
BOOST_REQUIRE(sats[phase].size() == s_opm[phase].size());
for (size_t i = 0; i < s_opm[phase].size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(sats[phase][i], s_opm[phase][i], 100.*reltol);
BOOST_CHECK_CLOSE(sats[phase][i], s_ecl[phase][i], reltol_ecl);
}
std::cout << std::endl;
}
const auto& rs = comp.rs();
const std::vector<double> rs_opm { // opm
74.624983020822540, 74.659590408801634, 74.694380353364522, 74.729353362649505,
74.764509945812975, 74.799850613032362, 74.835375875509555, 74.87108624547416,
74.906982236186707, 75.088917653469309, 52.5, 57.5,
62.5, 67.5, 72.5, 76.45954840804761,
76.70621044909619, 76.952877357524045, 77.199549133522638, 77.446225777283587};
const std::vector<double> rs_ecl { // eclipse
74.625114, 74.659706, 74.694481, 74.729439,
74.764580, 74.799904, 74.835419, 74.875252,
74.968628, 75.088951, 52.500000, 57.500000,
62.500000, 67.500000, 72.500000, 76.168388,
76.349953, 76.531532, 76.713142, 76.894775,};
const auto& rv = comp.rv();
const std::vector<double> rv_opm { // opm
2.50e-6, 7.50e-6, 1.25e-5, 1.75e-5,
2.25e-5, 2.75e-5, 3.25e-5, 3.75e-5,
4.25e-5, 2.51158386e-4, 2.52203372e-4, 5.75e-5,
6.25e-5, 6.75e-5, 7.25e-5, 7.75e-5,
8.25e-5, 8.75e-5, 9.25e-5, 9.75e-5};
const std::vector<double> rv_ecl { // eclipse
0.24999999E-05, 0.74999998E-05, 0.12500000E-04, 0.17500000E-04,
0.22500000E-04, 0.27500000E-04, 0.32500000E-04, 0.37500002E-04,
0.42500000E-04, 0.25115837E-03, 0.25220393E-03, 0.57500001E-04,
0.62500003E-04, 0.67499997E-04, 0.72499999E-04, 0.77500001E-04,
0.82500002E-04, 0.87499997E-04, 0.92499999E-04, 0.97500000E-04};
for (size_t i = 0; i < rv_opm.size(); ++i) {
//std::cout << std::setprecision(10) << sats[phase][i] << '\n';
BOOST_CHECK_CLOSE(rs[i], rs_opm[i], 100*reltol);
BOOST_CHECK_CLOSE(rs[i], rs_ecl[i], reltol_ecl);
BOOST_CHECK_CLOSE(rv[i], rv_opm[i], 100.*reltol);
BOOST_CHECK_CLOSE(rv[i], rv_ecl[i], reltol_ecl);
}
}
BOOST_AUTO_TEST_SUITE_END()

264
tests/test_parallel_linearsolver.cpp Normal file
View File

@ -0,0 +1,264 @@
/*
Copyright 2014 Dr. Markus Blatt - HPC-Simulation-Software & Services
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // to suppress our messages when throwing
#define BOOST_TEST_MODULE OPM-ParallelIterativeSolverTest
#include <boost/test/unit_test.hpp>
// MPI header
#if HAVE_MPI
#include <mpi.h>
#include <dune/common/parallel/indexset.hh>
#include <dune/common/parallel/communicator.hh>
#include <dune/common/parallel/remoteindices.hh>
#include <dune/common/mpicollectivecommunication.hh>
#include <dune/common/collectivecommunication.hh>
#include <dune/istl/owneroverlapcopy.hh>
#include <opm/core/linalg/ParallelIstlInformation.hpp>
#else
#error "This file needs to compiled with MPI support!"
#endif
#include <opm/core/linalg/LinearSolverFactory.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <memory>
#include <cstdlib>
#include <string>
struct MPIFixture {
MPIFixture()
{
int m_argc = boost::unit_test::framework::master_test_suite().argc;
char** m_argv = boost::unit_test::framework::master_test_suite().argv;
MPI_Init(&m_argc, &m_argv);
}
~MPIFixture()
{
MPI_Finalize();
}
};
BOOST_GLOBAL_FIXTURE(MPIFixture);
struct MyMatrix
{
MyMatrix(std::size_t rows, std::size_t nnz)
: data(nnz, 0.0), rowStart(rows+1, -1),
colIndex(nnz, -1)
{}
MyMatrix()
: data(), rowStart(), colIndex()
{}
std::vector<double> data;
std::vector<int> rowStart;
std::vector<int> colIndex;
};
typedef int LocalId;
typedef int GlobalId;
typedef Dune::OwnerOverlapCopyCommunication<GlobalId,LocalId> Communication;
typedef Dune::OwnerOverlapCopyAttributeSet GridAttributes;
typedef GridAttributes::AttributeSet GridFlag;
typedef Dune::ParallelLocalIndex<GridFlag> LocalIndex;
/// \brief Sets up a paralle Laplacian.
///
/// The process stores the unknowns with indices in the range [start, end).
/// As we use an overlapping domain decomposition, the process owns the indices
/// in the range [istart, iend]. If we would only used the indices in this range then
/// they form a partitioning of the whole index set.
/// \tparam I The type of the parallel index set (for convenience)
/// \param indexset The parallel index set for marking owner and copy region.
/// \param N The global number of unknowns of the system.
/// \param start The first index stored on this process
/// \param end One past the last index stored on this process
/// \param istart The first index that the process owns.
/// \param iend One past the last index the process owns.
template<class I>
std::shared_ptr<MyMatrix> create1DLaplacian(I& indexset, int N, int start, int end,
int istart, int iend)
{
indexset.beginResize();
MyMatrix* mm=new MyMatrix(end-start, (end-start)*3);
int nnz=0;
mm->rowStart[0]=0;
assert(start==0||start<istart);
assert(end==N||iend<end);
for(int row=start, localRow=0; row<end; row++, localRow++)
{
if(row<istart || row>=iend)
{
// We are in the overlap region of the grid
// therefore we setup the system such that
// right hand side will equal the left hand side
// of the linear system.
if(localRow>0)
{
mm->colIndex[nnz]=localRow-1;
mm->data[nnz++]=0;
}
mm->colIndex[nnz]=localRow;
mm->data[nnz++]=1.0;
indexset.add(row, LocalIndex(localRow, GridAttributes::copy, true));
if(localRow<end-1)
{
mm->colIndex[nnz]=localRow+1;
mm->data[nnz++]=0;
}
mm->rowStart[localRow+1]=nnz;
continue;
}
double dval=0;
if(row>0)
{
mm->colIndex[nnz]=localRow-1;
mm->data[nnz++]=-1;
dval+=1;
}
mm->colIndex[nnz]=localRow;
mm->data[nnz++]=2;//dval+(row<N-1);
if(row<N-1)
{
mm->colIndex[nnz]=localRow+1;
mm->data[nnz++]=-1;
dval+=1;
}
mm->rowStart[localRow+1]=nnz;
indexset.add(row, LocalIndex(localRow, GridAttributes::owner, true));
}
mm->data.resize(nnz);
mm->colIndex.resize(nnz);
indexset.endResize();
return std::shared_ptr<MyMatrix>(mm);
}
template<class O>
void createRandomVectors(O& pinfo, int NN, std::vector<double>& x, std::vector<double>& b,
const MyMatrix& mat)
{
x.resize(NN);
for(auto entry=x.begin(), end =x.end(); entry!=end; ++entry)
*entry=((double) (rand()%100))/10.0;
pinfo.copyOwnerToAll(x,x);
b.resize(NN);
// Construct the right hand side as b=A*x
std::fill(b.begin(), b.end(), 0.0);
for(std::size_t row=0; row<mat.rowStart.size()-1; ++row)
{
for(int i=mat.rowStart[row], end=mat.rowStart[row+1]; i!=end; ++i)
{
b[row]+= mat.data[i]*x[mat.colIndex[i]];
}
}
pinfo.copyOwnerToAll(b,b);
}
void run_test(const Opm::parameter::ParameterGroup& param)
{
int N=100;
int procs, rank;
MPI_Comm_rank(MPI_COMM_WORLD, &rank);
MPI_Comm_size(MPI_COMM_WORLD, &procs);
int n = N/procs; // number of unknowns per process
int bigger = N%procs; // number of process with n+1 unknows
int start, end, istart, iend;
// Compute owner region
if(rank<bigger) {
start = rank*(n+1);
end = start+(n+1);
}else{
start = bigger*(n+1) + (rank-bigger) * n;
end = start+n;
}
// Compute owner region
if(rank<bigger) {
istart = rank*(n+1);
iend = start+(n+1);
}else{
istart = bigger*(n+1) + (rank-bigger) * n;
iend = start+n;
}
// Compute overlap region
if(istart>0)
start = istart - 1;
else
start = istart;
if(iend<N)
end = iend + 1;
else
end = iend;
Opm::ParallelISTLInformation comm(MPI_COMM_WORLD);
auto mat = create1DLaplacian(*comm.indexSet(), N, start, end, istart, iend);
std::vector<double> x(end-start), b(end-start);
createRandomVectors(comm, end-start, x, b, *mat);
std::vector<double> exact(x);
std::fill(x.begin(), x.end(), 0.0);
Opm::LinearSolverFactory ls(param);
boost::any anyComm(comm);
ls.solve(b.size(), mat->data.size(), &(mat->rowStart[0]),
&(mat->colIndex[0]), &(mat->data[0]), &(b[0]),
&(x[0]), anyComm);
}
#ifdef HAVE_DUNE_ISTL
BOOST_AUTO_TEST_CASE(CGAMGTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("1"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(CGILUTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("0"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
run_test(param);
}
BOOST_AUTO_TEST_CASE(BiCGILUTest)
{
Opm::parameter::ParameterGroup param;
param.insertParameter(std::string("linsolver"), std::string("istl"));
param.insertParameter(std::string("linsolver_type"), std::string("2"));
param.insertParameter(std::string("linsolver_max_iterations"), std::string("200"));
run_test(param);
}
#endif

64
tests/test_regionmapping.cpp Normal file
View File

@ -0,0 +1,64 @@
/*
Copyright 2014 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"
/* --- Boost.Test boilerplate --- */
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE UnitsTest
#include <boost/test/unit_test.hpp>
#include <boost/test/floating_point_comparison.hpp>
/* --- our own headers --- */
#include <opm/core/utility/RegionMapping.hpp>
BOOST_AUTO_TEST_SUITE ()
BOOST_AUTO_TEST_CASE (RegionMapping)
{
// 0 1 2 3 4 5 6 7 8
std::vector<int> regions = { 2, 5, 2, 4, 2, 7, 6, 3, 6 };
Opm::RegionMapping<> rm(regions);
for (size_t i = 0; i < regions.size(); ++i) {
BOOST_CHECK_EQUAL(rm.region(i), regions[i]);
}
std::vector<int> region_ids = { 2, 3, 4, 5, 6, 7 };
std::vector< std::vector<int> > region_cells = { { 0, 2, 4 }, { 7 }, { 3 }, { 1 }, { 6, 8 }, { 5 } };
BOOST_REQUIRE_EQUAL(rm.numRegions(), region_ids.size());
for (size_t i = 0; i < region_ids.size(); ++i) {
auto cells = rm.cells(region_ids[i]);
BOOST_REQUIRE_EQUAL(std::distance(cells.begin(), cells.end()), region_cells[i].size());
size_t count = 0;
for (auto iter = cells.begin(); iter != cells.end(); ++iter, ++count) {
BOOST_CHECK_EQUAL(*iter, region_cells[i][count]);
}
}
}
BOOST_AUTO_TEST_SUITE_END()