OilfieldToolsNet/opm-simulators
4
0
Fork 0
mirror of https://github.com/OPM/opm-simulators.git synced 2025-02-25 18:55:30 -06:00
Code Issues Packages Projects Releases Wiki Activity

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

Browse Source 
Resolved Conflicts:
	opm/core/props/BlackoilPropertiesFromDeck.cpp
	opm/core/props/rock/RockFromDeck.hpp
	opm/core/props/satfunc/SaturationPropsFromDeck.hpp
	opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp
This commit is contained in:
Markus Blatt 2014-02-19 15:22:18 +01:00
commit 21864388bd
26 changed files with 2132 additions and 330 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

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

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

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

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

@ -27,6 +27,9 @@
#include <opm/core/props/satfunc/SaturationPropsFromDeck.hpp>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/utility/parameters/ParameterGroup.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <memory>
struct UnstructuredGrid;
@ -44,7 +47,15 @@ namespace Opm
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(const EclipseGridParser& deck,
BlackoilPropertiesFromDeck(const EclipseGridParser &deck,
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck and grid.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid, bool init_rock=true );
/// Initialize from deck, grid and parameters.
@ -63,6 +74,22 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock=true);
/// Initialize from deck, grid and parameters.
/// \param[in] deck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
/// \param[in] param Parameters. Accepted parameters include:
/// pvt_tab_size (200) number of uniform sample points for dead-oil pvt tables.
/// sat_tab_size (200) number of uniform sample points for saturation tables.
/// threephase_model("simple") three-phase relperm model (accepts "simple" and "stone2").
/// For both size parameters, a 0 or negative value indicates that no spline fitting is to
/// be done, and the input fluid data used directly for linear interpolation.
BlackoilPropertiesFromDeck(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const parameter::ParameterGroup& param,
bool init_rock=true);
/// Destructor.
virtual ~BlackoilPropertiesFromDeck();
@ -202,6 +229,23 @@ namespace Opm
const parameter::ParameterGroup& param,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
bool init_rock);
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cart_dims,
T begin_cell_centroids,
int dimension,
const parameter::ParameterGroup& param,
bool init_rock);
RockFromDeck rock_;
BlackoilPvtProperties pvt_;
std::unique_ptr<SaturationPropsInterface> satprops_;

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

@ -24,13 +24,13 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
namespace Opm
{
/// Looks at presence of WATER, OIL and GAS keywords in state object
/// to determine active phases.
inline PhaseUsage phaseUsageFromDeck(Opm::EclipseStateConstPtr eclipseState)
@ -105,6 +105,43 @@ namespace Opm
return pu;
}
/// Looks at presence of WATER, OIL and GAS keywords in deck
/// to determine active phases.
inline PhaseUsage phaseUsageFromDeck(Opm::DeckConstPtr newParserDeck)
{
PhaseUsage pu;
std::fill(pu.phase_used, pu.phase_used + BlackoilPhases::MaxNumPhases, 0);
// Discover phase usage.
if (newParserDeck->hasKeyword("WATER")) {
pu.phase_used[BlackoilPhases::Aqua] = 1;
}
if (newParserDeck->hasKeyword("OIL")) {
pu.phase_used[BlackoilPhases::Liquid] = 1;
}
if (newParserDeck->hasKeyword("GAS")) {
pu.phase_used[BlackoilPhases::Vapour] = 1;
}
pu.num_phases = 0;
for (int i = 0; i < BlackoilPhases::MaxNumPhases; ++i) {
pu.phase_pos[i] = pu.num_phases;
pu.num_phases += pu.phase_used[i];
}
// Only 2 or 3 phase systems handled.
if (pu.num_phases < 2 || pu.num_phases > 3) {
OPM_THROW(std::runtime_error, "Cannot handle cases with " << pu.num_phases << " phases.");
}
// We need oil systems, since we do not support the keywords needed for
// water-gas systems.
if (!pu.phase_used[BlackoilPhases::Liquid]) {
OPM_THROW(std::runtime_error, "Cannot handle cases with no OIL, i.e. water-gas systems.");
}
return pu;
}
}
#endif // OPM_PHASEUSAGEFROMDECK_HEADER_INCLUDED

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

@ -25,6 +25,9 @@
#include <opm/core/utility/ErrorMacros.hpp>
#include <opm/core/utility/linearInterpolation.hpp>
#include <opm/parser/eclipse/Utility/RocktabTable.hpp>
#include <opm/parser/eclipse/Utility/RockTable.hpp>
#include <iostream>
namespace Opm
@ -65,6 +68,44 @@ namespace Opm
}
}
RockCompressibility::RockCompressibility(Opm::DeckConstPtr newParserDeck)
: pref_(0.0),
rock_comp_(0.0)
{
if (newParserDeck->hasKeyword("ROCKTAB")) {
Opm::DeckKeywordConstPtr rtKeyword = newParserDeck->getKeyword("ROCKTAB");
if (rtKeyword->size() != 1)
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCKTAB.");
// the number of colums of the "ROCKTAB" keyword
// depends on the presence of the "RKTRMDIR"
// keyword. Messy stuff...
bool isDirectional = newParserDeck->hasKeyword("RKTRMDIR");
if (isDirectional)
{
// well, okay. we don't support non-isotropic
// transmissibility multipliers yet
OPM_THROW(std::runtime_error, "Support for non-isotropic "
"transmissibility multipliers is not implemented yet.");
};
Opm::RocktabTable rocktabTable(rtKeyword, isDirectional);
p_ = rocktabTable.getPressureColumn();
poromult_ = rocktabTable.getPoreVolumeMultiplierColumn();
transmult_ = rocktabTable.getTransmissibilityMultiplierColumn();
} else if (newParserDeck->hasKeyword("ROCK")) {
Opm::RockTable rockTable(newParserDeck->getKeyword("ROCK"));
if (rockTable.numRows() != 1)
OPM_THROW(std::runtime_error, "Can only handle a single region in ROCK.");
pref_ = rockTable.getPressureColumn()[0];
rock_comp_ = rockTable.getCompressibilityColumn()[0];
} else {
std::cout << "**** warning: no rock compressibility data found in deck (ROCK or ROCKTAB)." << std::endl;
}
}
bool RockCompressibility::isActive() const
{
return !p_.empty() || (rock_comp_ != 0.0);

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

@ -20,6 +20,8 @@
#ifndef OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#define OPM_ROCKCOMPRESSIBILITY_HEADER_INCLUDED
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
namespace Opm
@ -35,6 +37,10 @@ namespace Opm
/// Looks for the keywords ROCK and ROCKTAB.
RockCompressibility(const EclipseGridParser& deck);
/// Construct from input deck.
/// Looks for the keywords ROCK and ROCKTAB.
RockCompressibility(Opm::DeckConstPtr newParserDeck);
/// Construct from parameters.
/// Accepts the following parameters (with defaults).
/// rock_compressibility_pref (100.0) [given in bar]

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

@ -21,6 +21,9 @@
#include "config.h"
#include <opm/core/props/rock/RockFromDeck.hpp>
#include <opm/core/grid.h>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <array>
namespace Opm
@ -37,6 +40,10 @@ namespace Opm
PermeabilityKind fillTensor(const EclipseGridParser& parser,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap);
PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap);
} // anonymous namespace
@ -77,6 +84,17 @@ namespace Opm
assignPermeability(deck, number_of_cells, global_cell, cart_dims, perm_threshold);
}
void RockFromDeck::init(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell,
const int* cart_dims)
{
assignPorosity(newParserDeck, number_of_cells, global_cell);
permfield_valid_.assign(number_of_cells, false);
const double perm_threshold = 0.0; // Maybe turn into parameter?
assignPermeability(newParserDeck, number_of_cells, global_cell, cart_dims,
perm_threshold);
}
void RockFromDeck::assignPorosity(const EclipseGridParser& parser,
int number_of_cells, const int* global_cell)
@ -91,6 +109,20 @@ namespace Opm
}
}
void RockFromDeck::assignPorosity(Opm::DeckConstPtr newParserDeck,
int number_of_cells, const int* global_cell)
{
porosity_.assign(number_of_cells, 1.0);
if (newParserDeck->hasKeyword("PORO")) {
const std::vector<double>& poro = newParserDeck->getKeyword("PORO")->getSIDoubleData();
for (int c = 0; c < int(porosity_.size()); ++c) {
const int deck_pos = (global_cell == NULL) ? c : global_cell[c];
assert(0 <= c && c < (int) porosity_.size());
assert(0 <= deck_pos && deck_pos < (int) poro.size());
porosity_[c] = poro[deck_pos];
}
}
}
void RockFromDeck::assignPermeability(const EclipseGridParser& parser,
@ -147,6 +179,61 @@ namespace Opm
}
}
void RockFromDeck::assignPermeability(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const int* cartdims,
double perm_threshold)
{
const int dim = 3;
const int num_global_cells = cartdims[0]*cartdims[1]*cartdims[2];
const int nc = number_of_cells;
assert(num_global_cells > 0);
permeability_.assign(dim * dim * nc, 0.0);
std::vector<const std::vector<double>*> tensor;
tensor.reserve(10);
const std::vector<double> zero(num_global_cells, 0.0);
tensor.push_back(&zero);
std::array<int,9> kmap;
PermeabilityKind pkind = fillTensor(newParserDeck, tensor, kmap);
if (pkind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid permeability field.");
}
// Assign permeability values only if such values are
// given in the input deck represented by 'newParserDeck'. In
// other words: Don't set any (arbitrary) default values.
// It is infinitely better to experience a reproducible
// crash than subtle errors resulting from a (poorly
// chosen) default value...
//
if (tensor.size() > 1) {
const int* gc = global_cell;
int off = 0;
for (int c = 0; c < nc; ++c, off += dim*dim) {
// SharedPermTensor K(dim, dim, &permeability_[off]);
int kix = 0;
const int glob = (gc == NULL) ? c : gc[c];
for (int i = 0; i < dim; ++i) {
for (int j = 0; j < dim; ++j, ++kix) {
// K(i,j) = (*tensor[kmap[kix]])[glob];
permeability_[off + kix] = (*tensor[kmap[kix]])[glob];
}
// K(i,i) = std::max(K(i,i), perm_threshold);
permeability_[off + 3*i + i] = std::max(permeability_[off + 3*i + i], perm_threshold);
}
permfield_valid_[c] = std::vector<unsigned char>::value_type(1);
}
}
}
namespace {
@ -216,6 +303,72 @@ namespace Opm
return retval;
}
/// @brief
/// Classify and verify a given permeability specification
/// from a structural point of view. In particular, we
/// verify that there are no off-diagonal permeability
/// components such as @f$k_{xy}@f$ unless the
/// corresponding diagonal components are known as well.
///
/// @param newParserDeck [in]
/// An Eclipse data parser capable of answering which
/// permeability components are present in a given input
/// deck.
///
/// @return
/// An enum value with the following possible values:
/// ScalarPerm only one component was given.
/// DiagonalPerm more than one component given.
/// TensorPerm at least one cross-component given.
/// None no components given.
/// Invalid invalid set of components given.
PermeabilityKind classifyPermeability(Opm::DeckConstPtr newParserDeck)
{
const bool xx = newParserDeck->hasKeyword("PERMX" );
const bool xy = newParserDeck->hasKeyword("PERMXY");
const bool xz = newParserDeck->hasKeyword("PERMXZ");
const bool yx = newParserDeck->hasKeyword("PERMYX");
const bool yy = newParserDeck->hasKeyword("PERMY" );
const bool yz = newParserDeck->hasKeyword("PERMYZ");
const bool zx = newParserDeck->hasKeyword("PERMZX");
const bool zy = newParserDeck->hasKeyword("PERMZY");
const bool zz = newParserDeck->hasKeyword("PERMZ" );
int num_cross_comp = xy + xz + yx + yz + zx + zy;
int num_comp = xx + yy + zz + num_cross_comp;
PermeabilityKind retval = None;
if (num_cross_comp > 0) {
retval = TensorPerm;
} else {
if (num_comp == 1) {
retval = ScalarPerm;
} else if (num_comp >= 2) {
retval = DiagonalPerm;
}
}
bool ok = true;
if (num_comp > 0) {
// At least one tensor component specified on input.
// Verify that any remaining components are OK from a
// structural point of view. In particular, there
// must not be any cross-components (e.g., k_{xy})
// unless the corresponding diagonal component (e.g.,
// k_{xx}) is present as well...
//
ok = xx || !(xy || xz || yx || zx) ;
ok = ok && (yy || !(yx || yz || xy || zy));
ok = ok && (zz || !(zx || zy || xz || yz));
}
if (!ok) {
retval = Invalid;
}
return retval;
}
/// @brief
/// Copy isotropic (scalar) permeability to other diagonal
@ -345,6 +498,106 @@ namespace Opm
return kind;
}
/// @brief
/// Extract pointers to appropriate tensor components from
/// input deck. The permeability tensor is, generally,
/// @code
/// [ kxx kxy kxz ]
/// K = [ kyx kyy kyz ]
/// [ kzx kzy kzz ]
/// @endcode
/// We store these values in a linear array using natural
/// ordering with the column index cycling the most rapidly.
/// In particular we use the representation
/// @code
/// [ 0 1 2 3 4 5 6 7 8 ]
/// K = [ kxx, kxy, kxz, kyx, kyy, kyz, kzx, kzy, kzz ]
/// @endcode
/// Moreover, we explicitly enforce symmetric tensors by
/// assigning
/// @code
/// 3 1 6 2 7 5
/// kyx = kxy, kzx = kxz, kzy = kyz
/// @endcode
/// However, we make no attempt at enforcing positive
/// definite tensors.
///
/// @param [in] parser
/// An Eclipse data parser capable of answering which
/// permeability components are present in a given input
/// deck as well as retrieving the numerical value of each
/// permeability component in each grid cell.
///
/// @param [out] tensor
/// @param [out] kmap
PermeabilityKind fillTensor(Opm::DeckConstPtr newParserDeck,
std::vector<const std::vector<double>*>& tensor,
std::array<int,9>& kmap)
{
PermeabilityKind kind = classifyPermeability(newParserDeck);
if (kind == Invalid) {
OPM_THROW(std::runtime_error, "Invalid set of permeability fields given.");
}
assert(tensor.size() == 1);
for (int i = 0; i < 9; ++i) { kmap[i] = 0; }
enum { xx, xy, xz, // 0, 1, 2
yx, yy, yz, // 3, 4, 5
zx, zy, zz }; // 6, 7, 8
// -----------------------------------------------------------
// 1st row: [kxx, kxy, kxz]
if (newParserDeck->hasKeyword("PERMX" )) {
kmap[xx] = tensor.size();
tensor.push_back(&newParserDeck->getKeyword("PERMX")->getSIDoubleData());
setScalarPermIfNeeded(kmap, xx, yy, zz);
}
if (newParserDeck->hasKeyword("PERMXY")) {
kmap[xy] = kmap[yx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMXY")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMXZ")) {
kmap[xz] = kmap[zx] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMXZ")->getSIDoubleData());
}
// -----------------------------------------------------------
// 2nd row: [kyx, kyy, kyz]
if (newParserDeck->hasKeyword("PERMYX")) {
kmap[yx] = kmap[xy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMYX")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMY" )) {
kmap[yy] = tensor.size();
tensor.push_back(&newParserDeck->getKeyword("PERMY")->getSIDoubleData());
setScalarPermIfNeeded(kmap, yy, zz, xx);
}
if (newParserDeck->hasKeyword("PERMYZ")) {
kmap[yz] = kmap[zy] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMYZ")->getSIDoubleData());
}
// -----------------------------------------------------------
// 3rd row: [kzx, kzy, kzz]
if (newParserDeck->hasKeyword("PERMZX")) {
kmap[zx] = kmap[xz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMZX")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMZY")) {
kmap[zy] = kmap[yz] = tensor.size(); // Enforce symmetry.
tensor.push_back(&newParserDeck->getKeyword("PERMZY")->getSIDoubleData());
}
if (newParserDeck->hasKeyword("PERMZ" )) {
kmap[zz] = tensor.size();
tensor.push_back(&newParserDeck->getKeyword("PERMZ")->getSIDoubleData());
setScalarPermIfNeeded(kmap, zz, xx, yy);
}
return kind;
}
} // anonymous namespace
} // namespace Opm

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

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

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

@ -27,6 +27,9 @@
#include <opm/core/props/satfunc/SatFuncStone2.hpp>
#include <opm/core/props/satfunc/SatFuncSimple.hpp>
#include <opm/core/props/satfunc/SatFuncGwseg.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <vector>
struct UnstructuredGrid;
@ -82,6 +85,38 @@ namespace Opm
int dimensions,
const int samples);
/// Initialize from deck and grid.
/// \param[in] newParserDeck Deck input parser
/// \param[in] grid Grid to which property object applies, needed for the
/// mapping from cell indices (typically from a processed grid)
/// to logical cartesian indices consistent with the deck.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
void init(Opm::DeckConstPtr newParserDeck,
const UnstructuredGrid& grid,
const int samples);
/// Initialize from deck and grid.
/// \param[in] newParserDeck Deck input parser
/// \param[in] number_of_cells The number of cells of the grid to which property
/// object applies, needed for the
/// mapping from cell indices (typically from a processed
/// grid) to logical cartesian indices consistent with the
/// deck.
/// \param[in] global_cell The mapping from local cell indices of the grid to
/// global cell indices used in the deck.
/// \param[in] begin_cell_centroids Pointer to the first cell_centroid of the grid.
/// \param[in] dimensions The dimensions of the grid.
/// \param[in] samples Number of uniform sample points for saturation tables.
/// NOTE: samples will only be used with the SatFuncSetUniform template argument.
template<class T>
void init(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const int samples);
/// \return P, the number of phases.
int numPhases() const;
@ -123,37 +158,92 @@ namespace Opm
const int* cells,
double* smin,
double* smax) const;
/// Update saturation state for the hysteresis tracking
/// \param[in] n Number of data points.
/// \param[in] s Array of nP saturation values.
void updateSatHyst(const int n,
const int* cells,
const double* s);
private:
PhaseUsage phase_usage_;
std::vector<SatFuncSet> satfuncset_;
std::vector<int> cell_to_func_; // = SATNUM - 1
struct { // End point scaling parameters
std::vector<double> swl_;
std::vector<double> swcr_;
std::vector<double> swu_;
std::vector<double> sowcr_;
std::vector<double> krw_;
std::vector<double> krwr_;
std::vector<double> kro_;
std::vector<double> krorw_;
} eps_;
bool do_eps_; // ENDSCALE is active
bool do_3pt_; // SCALECRS: YES~true NO~false
bool do_hyst_; // Keywords ISWL etc detected
std::vector<EPSTransforms> eps_transf_;
std::vector<EPSTransforms> eps_transf_hyst_;
std::vector<SatHyst> sat_hyst_;
typedef SatFuncSet Funcs;
const Funcs& funcForCell(const int cell) const;
template<class T>
void initEPS(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
void relpermEPS(const double *s, const int cell, double *kr, double *dkrds= 0) const;
int dimensions);
template<class T>
void initEPSHyst(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(const EclipseGridParser& deck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
template<class T>
void initEPS(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSHyst(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions);
template<class T>
void initEPSKey(Opm::DeckConstPtr newParserDeck,
int number_of_cells,
const int* global_cell,
const T& begin_cell_centroids,
int dimensions,
const std::string& keyword,
std::vector<double>& scaleparam);
void initEPSParam(const int cell,
EPSTransforms::Transform& data,
const bool oil,
const double sl_tab,
const double scr_tab,
const double su_tab,
const double sxcr_tab,
const double s0_tab,
const double krsr_tab,
const double krmax_tab,
const std::vector<double>& sl,
const std::vector<double>& scr,
const std::vector<double>& su,
const std::vector<double>& sxcr,
const std::vector<double>& s0,
const std::vector<double>& krsr,
const std::vector<double>& krmax);
bool columnIsMasked_(Opm::DeckConstPtr newParserDeck,
const std::string &keywordName,
int columnIdx)
{ return newParserDeck->getKeyword(keywordName)->getRecord(0)->getItem(0)->getSIDouble(0) != -1.0; }
};

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

File diff suppressed because it is too large Load Diff

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

@ -33,6 +33,8 @@
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
#include <opm/parser/eclipse/Utility/EquilWrapper.hpp>
#include <iostream>
#include <cmath>
@ -547,6 +549,103 @@ namespace Opm
}
/// Initialize a state from input deck.
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
Opm::DeckConstPtr newParserDeck,
const double gravity,
State& state)
{
const int num_phases = props.numPhases();
const PhaseUsage pu = phaseUsageFromDeck(newParserDeck);
if (num_phases != pu.num_phases) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): user specified property object with " << num_phases << " phases, "
"found " << pu.num_phases << " phases in deck.");
}
state.init(grid, num_phases);
if (newParserDeck->hasKeyword("EQUIL")) {
if (num_phases != 2) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
}
if (pu.phase_used[BlackoilPhases::Vapour]) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas).");
}
// Set saturations depending on oil-water contact.
EquilWrapper equil(newParserDeck->getKeyword("EQUIL"));
if (equil.numRegions() != 1) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): No region support yet.");
}
const double woc = equil.waterOilContactDepth(0);
initWaterOilContact(grid, props, woc, WaterBelow, state);
// Set pressure depending on densities and depths.
const double datum_z = equil.datumDepth(0);
const double datum_p = equil.datumDepthPressure(0);
initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state);
} else if (newParserDeck->hasKeyword("PRESSURE")) {
// Set saturations from SWAT/SGAS, pressure from PRESSURE.
std::vector<double>& s = state.saturation();
std::vector<double>& p = state.pressure();
const std::vector<double>& p_deck = newParserDeck->getKeyword("PRESSURE")->getSIDoubleData();
const int num_cells = grid.number_of_cells;
if (num_phases == 2) {
if (!pu.phase_used[BlackoilPhases::Aqua]) {
// oil-gas: we require SGAS
if (!newParserDeck->hasKeyword("SGAS")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS keyword in 2-phase init");
}
const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c + gpos] = sg_deck[c_deck];
s[2*c + opos] = 1.0 - sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
// water-oil or water-gas: we require SWAT
if (!newParserDeck->hasKeyword("SWAT")) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SWAT keyword in 2-phase init");
}
const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int nwpos = (wpos + 1) % 2;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[2*c + wpos] = sw_deck[c_deck];
s[2*c + nwpos] = 1.0 - sw_deck[c_deck];
p[c] = p_deck[c_deck];
}
}
} else if (num_phases == 3) {
const bool has_swat_sgas = newParserDeck->hasKeyword("SWAT") && newParserDeck->hasKeyword("SGAS");
if (!has_swat_sgas) {
OPM_THROW(std::runtime_error, "initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init.");
}
const int wpos = pu.phase_pos[BlackoilPhases::Aqua];
const int gpos = pu.phase_pos[BlackoilPhases::Vapour];
const int opos = pu.phase_pos[BlackoilPhases::Liquid];
const std::vector<double>& sw_deck = newParserDeck->getKeyword("SWAT")->getSIDoubleData();
const std::vector<double>& sg_deck = newParserDeck->getKeyword("SGAS")->getSIDoubleData();
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
s[3*c + wpos] = sw_deck[c_deck];
s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
s[3*c + gpos] = sg_deck[c_deck];
p[c] = p_deck[c_deck];
}
} else {
OPM_THROW(std::runtime_error, "initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases.");
}
} else {
OPM_THROW(std::runtime_error, "initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS.");
}
// Finally, init face pressures.
initFacePressure(grid, state);
}
/// Initialize a state from input deck.
template <class Props, class State>
void initStateFromDeck(const UnstructuredGrid& grid,
@ -666,10 +765,6 @@ namespace Opm
begin_cell_centroids, state);
}
/// Initialize surface volume from pressure and saturation by z = As.
/// Here saturation is used as an intial guess for z in the
/// computation of A.
@ -904,6 +999,39 @@ namespace Opm
}
}
/// Initialize a blackoil state from input deck.
template <class Props, class State>
void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
const Props& props,
Opm::DeckConstPtr newParserDeck,
const double gravity,
State& state)
{
initStateFromDeck(grid, props, newParserDeck, gravity, state);
if (newParserDeck->hasKeyword("RS")) {
const std::vector<double>& rs_deck = newParserDeck->getKeyword("RS")->getSIDoubleData();
const int num_cells = grid.number_of_cells;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
state.gasoilratio()[c] = rs_deck[c_deck];
}
initBlackoilSurfvolUsingRSorRV(grid, props, state);
computeSaturation(props,state);
} else if (newParserDeck->hasKeyword("RV")){
const std::vector<double>& rv_deck = newParserDeck->getKeyword("RV")->getSIDoubleData();
const int num_cells = grid.number_of_cells;
for (int c = 0; c < num_cells; ++c) {
int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
state.rv()[c] = rv_deck[c_deck];
}
initBlackoilSurfvolUsingRSorRV(grid, props, state);
computeSaturation(props,state);
}
else {
OPM_THROW(std::runtime_error, "Temporarily, we require the RS or the RV field.");
}
}
} // namespace Opm

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

@ -19,10 +19,66 @@
#include "config.h"
#include <opm/core/wells/WellCollection.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <boost/lexical_cast.hpp>
#include <memory>
namespace Opm
{
void WellCollection::addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage) {
WellsGroupInterface* fieldNode = findNode(fieldGroup->name());
if (fieldNode) {
OPM_THROW(std::runtime_error, "Trying to add FIELD node, but this already exists. Can only have one FIELD node.");
}
roots_.push_back(createGroupWellsGroup(fieldGroup, timeStep, phaseUsage));
}
void WellCollection::addGroup(GroupConstPtr groupChild, std::string parent_name,
size_t timeStep, const PhaseUsage& phaseUsage) {
WellsGroupInterface* parent = findNode(parent_name);
if (!parent) {
OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent_name << ", but this does not exist in the WellCollection.");
}
if (findNode(groupChild->name())) {
OPM_THROW(std::runtime_error, "Trying to add child group named " << groupChild->name() << ", but this group is already in the WellCollection.");
}
std::shared_ptr<WellsGroupInterface> child = createGroupWellsGroup(groupChild, timeStep, phaseUsage);
WellsGroup* parent_as_group = static_cast<WellsGroup*> (parent);
if (!parent_as_group) {
OPM_THROW(std::runtime_error, "Trying to add child group to group named " << parent->name() << ", but it's not a group.");
}
parent_as_group->addChild(child);
child->setParent(parent);
}
void WellCollection::addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage) {
WellsGroupInterface* parent = findNode(wellChild->getGroupName(timeStep));
if (!parent) {
OPM_THROW(std::runtime_error, "Trying to add well " << wellChild->name() << " Step: " << boost::lexical_cast<std::string>(timeStep) << " to group named " << wellChild->getGroupName(timeStep) << ", but this group does not exist in the WellCollection.");
}
std::shared_ptr<WellsGroupInterface> child = createWellWellsGroup(wellChild, timeStep, phaseUsage);
WellsGroup* parent_as_group = static_cast<WellsGroup*> (parent);
if (!parent_as_group) {
OPM_THROW(std::runtime_error, "Trying to add well to group named " << wellChild->getGroupName(timeStep) << ", but it's not a group.");
}
parent_as_group->addChild(child);
leaf_nodes_.push_back(static_cast<WellNode*>(child.get()));
child->setParent(parent);
}
void WellCollection::addChild(const std::string& child_name,
const std::string& parent_name,
@ -33,6 +89,7 @@ namespace Opm
roots_.push_back(createWellsGroup(parent_name, deck));
parent = roots_[roots_.size() - 1].get();
}
std::shared_ptr<WellsGroupInterface> child;
for (size_t i = 0; i < roots_.size(); ++i) {
@ -47,6 +104,7 @@ namespace Opm
break;
}
}
if (!child.get()) {
child = createWellsGroup(child_name, deck);
}
@ -65,7 +123,6 @@ namespace Opm
}
const std::vector<WellNode*>& WellCollection::getLeafNodes() const {
return leaf_nodes_;
}

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

@ -23,10 +23,15 @@
#define OPM_WELLCOLLECTION_HPP
#include <vector>
#include <memory>
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/core/grid.h>
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <memory>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
namespace Opm
{
@ -34,6 +39,14 @@ namespace Opm
class WellCollection
{
public:
void addField(GroupConstPtr fieldGroup, size_t timeStep, const PhaseUsage& phaseUsage);
void addWell(WellConstPtr wellChild, size_t timeStep, const PhaseUsage& phaseUsage);
void addGroup(GroupConstPtr groupChild, std::string parent_name,
size_t timeStep, const PhaseUsage& phaseUsage);
/// Adds and creates if necessary the child to the collection
/// and appends it to parent's children. Also adds and creates the parent
/// if necessary.

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

@ -74,7 +74,7 @@ namespace Opm
{
return parent_;
}
const std::string& WellsGroupInterface::name()
const std::string& WellsGroupInterface::name() const
{
return name_;
}
@ -1141,4 +1141,54 @@ namespace Opm
return return_value;
}
std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep, const PhaseUsage& phase_usage )
{
InjectionSpecification injection_specification;
ProductionSpecification production_specification;
if (group->isInjectionGroup(timeStep)) {
injection_specification.injector_type_ = toInjectorType(Phase::PhaseEnum2String(group->getInjectionPhase(timeStep)));
injection_specification.control_mode_ = toInjectionControlMode(GroupInjection::ControlEnum2String(group->getInjectionControlMode(timeStep)));
injection_specification.surface_flow_max_rate_ = group->getSurfaceMaxRate(timeStep);
injection_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep);
injection_specification.reinjection_fraction_target_ = group->getTargetReinjectFraction(timeStep);
injection_specification.voidage_replacment_fraction_ = group->getTargetVoidReplacementFraction(timeStep);
}
else if (group->isProductionGroup(timeStep)) {
production_specification.oil_max_rate_ = group->getOilTargetRate(timeStep);
production_specification.control_mode_ = toProductionControlMode(GroupProduction::ControlEnum2String(group->getProductionControlMode(timeStep)));
production_specification.water_max_rate_ = group->getWaterTargetRate(timeStep);
production_specification.gas_max_rate_ = group->getGasTargetRate(timeStep);
production_specification.liquid_max_rate_ = group->getLiquidTargetRate(timeStep);
production_specification.procedure_ = toProductionProcedure(GroupProductionExceedLimit::ActionEnum2String(group->getProductionExceedLimitAction(timeStep)));
production_specification.reservoir_flow_max_rate_ = group->getReservoirMaxRate(timeStep);
}
std::shared_ptr<WellsGroupInterface> wells_group(new WellsGroup(group->name(), production_specification, injection_specification, phase_usage));
return wells_group;
}
std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep, const PhaseUsage& phase_usage )
{
InjectionSpecification injection_specification;
ProductionSpecification production_specification;
if (well->isInjector(timeStep)) {
injection_specification.BHP_limit_ = well->getBHPLimit(timeStep);
injection_specification.injector_type_ = toInjectorType(WellInjector::Type2String(well->getInjectorType(timeStep)));
injection_specification.control_mode_ = toInjectionControlMode(WellInjector::ControlMode2String(well->getInjectorControlMode(timeStep)));
injection_specification.surface_flow_max_rate_ = well->getSurfaceInjectionRate(timeStep);
injection_specification.reservoir_flow_max_rate_ = well->getReservoirInjectionRate(timeStep);
production_specification.guide_rate_ = 0.0; // We know we're not a producer
}
else if (well->isProducer(timeStep)) {
production_specification.BHP_limit_ = well->getBHPLimit(timeStep);
production_specification.reservoir_flow_max_rate_ = well->getResVRate(timeStep);
production_specification.oil_max_rate_ = well->getOilRate(timeStep);
production_specification.control_mode_ = toProductionControlMode(WellProducer::ControlMode2String(well->getProducerControlMode(timeStep)));
production_specification.water_max_rate_ = well->getWaterRate(timeStep);
injection_specification.guide_rate_ = 0.0; // we know we're not an injector
}
std::shared_ptr<WellsGroupInterface> wells_group(new WellNode(well->name(), production_specification, injection_specification, phase_usage));
return wells_group;
}
}

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

@ -25,6 +25,10 @@
#include <opm/core/io/eclipse/EclipseGridParser.hpp>
#include <opm/core/grid.h>
#include <opm/core/props/BlackoilPhases.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <string>
#include <memory>
@ -58,7 +62,7 @@ namespace Opm
virtual ~WellsGroupInterface();
/// The unique identifier for the well or well group.
const std::string& name();
const std::string& name() const;
/// Production specifications for the well or well group.
const ProductionSpecification& prodSpec() const;
@ -403,9 +407,21 @@ namespace Opm
/// \param[in] name the name of the wells group.
/// \param[in] deck the deck from which to fetch information.
std::shared_ptr<WellsGroupInterface> createWellsGroup(const std::string& name,
const EclipseGridParser& deck);
const EclipseGridParser& deck);
/// Creates the WellsGroupInterface for the given well
/// \param[in] well the Well to construct object for
/// \param[in] timeStep the time step in question
/// \param[in] the phase usage
std::shared_ptr<WellsGroupInterface> createWellWellsGroup(WellConstPtr well, size_t timeStep,
const PhaseUsage& phase_usage );
/// Creates the WellsGroupInterface for the given Group
/// \param[in] group the Group to construct object for
/// \param[in] timeStep the time step in question
/// \param[in] the phase usage
std::shared_ptr<WellsGroupInterface> createGroupWellsGroup(GroupConstPtr group, size_t timeStep,
const PhaseUsage& phase_usage );
}
#endif /* OPM_WELLSGROUP_HPP */

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

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

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

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

3
tests/testBlackoilState1.DATA
View File

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

3
tests/testBlackoilState2.DATA
View File

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

91
tests/test_wellcollection.cpp Normal file
View File

@ -0,0 +1,91 @@
/*
Copyright 2012 SINTEF ICT, Applied Mathematics.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE WellCollectionTest
#include <boost/test/unit_test.hpp>
#include <opm/core/wells/WellCollection.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/Deck/Deck.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTreeNode.hpp>
using namespace Opm;
BOOST_AUTO_TEST_CASE(AddWellsAndGroupToCollection) {
ParserPtr parser(new Parser());
boost::filesystem::path scheduleFile("wells_group.data");
DeckConstPtr deck = parser->parseFile(scheduleFile.string());
EclipseStateConstPtr eclipseState(new EclipseState(deck));
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
GroupTreeNodePtr field=eclipseState->getSchedule()->getGroupTree(2)->getNode("FIELD");
GroupTreeNodePtr g1=eclipseState->getSchedule()->getGroupTree(2)->getNode("G1");
GroupTreeNodePtr g2=eclipseState->getSchedule()->getGroupTree(2)->getNode("G2");
WellCollection collection;
// Add groups to WellCollection
GroupConstPtr fieldGroup = eclipseState->getSchedule()->getGroup(field->name());
collection.addField(fieldGroup, 2, pu);
for (auto iter = field->begin(); iter != field->end(); ++iter) {
GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
collection.addGroup(childGroupNode, fieldGroup->name(), 2, pu);
}
GroupConstPtr g1Group = eclipseState->getSchedule()->getGroup(g1->name());
for (auto iter = g1->begin(); iter != g1->end(); ++iter) {
GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
collection.addGroup(childGroupNode, g1Group->name(), 2, pu);
}
GroupConstPtr g2Group = eclipseState->getSchedule()->getGroup(g2->name());
for (auto iter = g2->begin(); iter != g2->end(); ++iter) {
GroupConstPtr childGroupNode = eclipseState->getSchedule()->getGroup((*iter).second->name());
collection.addGroup(childGroupNode, g2Group->name(), 2, pu);
}
BOOST_CHECK_EQUAL("FIELD", collection.findNode("FIELD")->name());
BOOST_CHECK_EQUAL("FIELD", collection.findNode("G1")->getParent()->name());
BOOST_CHECK_EQUAL("FIELD", collection.findNode("G2")->getParent()->name());
// Add wells to WellCollection
WellCollection wellCollection;
std::vector<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
for (size_t i=0; i<wells.size(); i++) {
collection.addWell(wells[i], 2, pu);
}
BOOST_CHECK_EQUAL("G1", collection.findNode("INJ1")->getParent()->name());
BOOST_CHECK_EQUAL("G1", collection.findNode("INJ2")->getParent()->name());
BOOST_CHECK_EQUAL("G2", collection.findNode("PROD1")->getParent()->name());
BOOST_CHECK_EQUAL("G2", collection.findNode("PROD2")->getParent()->name());
}

108
tests/test_wellsgroup.cpp Normal file
View File

@ -0,0 +1,108 @@
/*
Copyright 2014 Statoil.
This file is part of the Open Porous Media project (OPM).
OPM is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
OPM is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <config.h>
#if HAVE_DYNAMIC_BOOST_TEST
#define BOOST_TEST_DYN_LINK
#endif
#define NVERBOSE // Suppress own messages when throw()ing
#define BOOST_TEST_MODULE WellsGroupTest
#include <memory>
#include <vector>
#include <boost/test/unit_test.hpp>
#include <boost/date_time/posix_time/posix_time.hpp>
#include <opm/core/props/phaseUsageFromDeck.hpp>
#include <opm/core/wells/WellsGroup.hpp>
#include <opm/parser/eclipse/Parser/Parser.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Group.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/GroupTreeNode.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well.hpp>
#include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
using namespace Opm;
BOOST_AUTO_TEST_CASE(ConstructGroupFromWell) {
ParserPtr parser(new Parser());
boost::filesystem::path scheduleFile("wells_group.data");
DeckConstPtr deck = parser->parseFile(scheduleFile.string());
EclipseStateConstPtr eclipseState(new EclipseState(deck));
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
std::vector<WellConstPtr> wells = eclipseState->getSchedule()->getWells();
for (size_t i=0; i<wells.size(); i++) {
WellConstPtr well = wells[i];
std::shared_ptr<WellsGroupInterface> wellsGroup = createWellWellsGroup(well, 2, pu);
BOOST_CHECK_EQUAL(well->name(), wellsGroup->name());
if (well->isInjector(2)) {
BOOST_CHECK_EQUAL(well->getSurfaceInjectionRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->injSpec().BHP_limit_);
BOOST_CHECK_EQUAL(well->getReservoirInjectionRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->prodSpec().guide_rate_);
}
if (well->isProducer(2)) {
BOOST_CHECK_EQUAL(well->getResVRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(well->getBHPLimit(2), wellsGroup->prodSpec().BHP_limit_);
BOOST_CHECK_EQUAL(well->getOilRate(2), wellsGroup->prodSpec().oil_max_rate_);
BOOST_CHECK_EQUAL(well->getWaterRate(2), wellsGroup->prodSpec().water_max_rate_);
BOOST_CHECK_EQUAL(0.0, wellsGroup->injSpec().guide_rate_);
}
}
}
BOOST_AUTO_TEST_CASE(ConstructGroupFromGroup) {
ParserPtr parser(new Parser());
boost::filesystem::path scheduleFile("wells_group.data");
DeckConstPtr deck = parser->parseFile(scheduleFile.string());
EclipseStateConstPtr eclipseState(new EclipseState(deck));
PhaseUsage pu = phaseUsageFromDeck(eclipseState);
std::vector<GroupTreeNodeConstPtr> nodes = eclipseState->getSchedule()->getGroupTree(2)->getNodes();
for (size_t i=0; i<nodes.size(); i++) {
GroupConstPtr group = eclipseState->getSchedule()->getGroup(nodes[i]->name());
std::shared_ptr<WellsGroupInterface> wellsGroup = createGroupWellsGroup(group, 2, pu);
BOOST_CHECK_EQUAL(group->name(), wellsGroup->name());
if (group->isInjectionGroup(2)) {
BOOST_CHECK_EQUAL(group->getSurfaceMaxRate(2), wellsGroup->injSpec().surface_flow_max_rate_);
BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->injSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(group->getTargetReinjectFraction(2), wellsGroup->injSpec().reinjection_fraction_target_);
BOOST_CHECK_EQUAL(group->getTargetVoidReplacementFraction(2), wellsGroup->injSpec().voidage_replacment_fraction_);
}
if (group->isProductionGroup(2)) {
BOOST_CHECK_EQUAL(group->getReservoirMaxRate(2), wellsGroup->prodSpec().reservoir_flow_max_rate_);
BOOST_CHECK_EQUAL(group->getGasTargetRate(2), wellsGroup->prodSpec().gas_max_rate_);
BOOST_CHECK_EQUAL(group->getOilTargetRate(2), wellsGroup->prodSpec().oil_max_rate_);
BOOST_CHECK_EQUAL(group->getWaterTargetRate(2), wellsGroup->prodSpec().water_max_rate_);
BOOST_CHECK_EQUAL(group->getLiquidTargetRate(2), wellsGroup->prodSpec().liquid_max_rate_);
}
}
}

31
tests/test_wellsmanager.cpp
View File

@ -207,7 +207,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(eclipseState, 0, Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
std::cout << "Checking new well structure, epoch 0" << std::endl;
@ -223,7 +223,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works) {
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(eclipseState, 1,Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
std::cout << "Checking new well structure, epoch 1" << std::endl;
@ -249,7 +249,7 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
Deck.setCurrentEpoch(0);
{
Opm::WellsManager wellsManager(eclipseState, 0, Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 0, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal(wellsManager.c_wells(), oldWellsManager.c_wells(),false));
@ -257,7 +257,15 @@ BOOST_AUTO_TEST_CASE(New_Constructor_Works_ExpandedData) {
Deck.setCurrentEpoch(1);
{
Opm::WellsManager wellsManager(eclipseState, 1,Deck, *gridManager.c_grid(), NULL);
Opm::WellsManager wellsManager(eclipseState, 1, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
}
Deck.setCurrentEpoch(2);
{
Opm::WellsManager wellsManager(eclipseState, 2, *gridManager.c_grid(), NULL);
Opm::WellsManager oldWellsManager(Deck, *gridManager.c_grid(), NULL);
BOOST_CHECK(wells_equal( wellsManager.c_wells(), oldWellsManager.c_wells(), true));
@ -303,3 +311,18 @@ BOOST_AUTO_TEST_CASE(ControlsEqual) {
}
BOOST_AUTO_TEST_CASE(WellHasSTOP_ExceptionIsThrown) {
Opm::EclipseGridParser Deck("wells_manager_data_wellSTOP.data");
Opm::GridManager gridManager(Deck);
Opm::ParserPtr parser(new Opm::Parser());
Opm::EclipseStateConstPtr eclipseState(new Opm::EclipseState(parser->parseFile("wells_manager_data_wellSTOP.data")));
Deck.setCurrentEpoch(0);
BOOST_CHECK_THROW( new Opm::WellsManager(eclipseState, 0, *gridManager.c_grid(), NULL), std::runtime_error );
}

68
tests/wells_group.data Executable file
View File

@ -0,0 +1,68 @@
OIL
GAS
WATER
DIMENS
10 10 5 /
GRID
DXV
10*1000.0 /
DYV
10*1000.0 /
DZV
10.0 20.0 30.0 10.0 5.0 /
DEPTHZ
121*2000
/
SCHEDULE
GRUPTREE
'G1' 'FIELD' /
'G2' 'FIELD' /
/
WELSPECS
'INJ1' 'G1' 1 1 8335 'GAS' /
'PROD1' 'G2' 10 10 8400 'OIL' /
/
TSTEP
14.0 /
/
WELSPECS
'INJ2' 'G1' 1 1 8335 'GAS' /
'PROD2' 'G2' 10 10 8400 'OIL' /
/
GCONINJE
'G1' GAS RATE 30000 /
/
GCONPROD
'G2' ORAT 10000 /
/
WCONINJE
'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 /
'INJ2' 'WATER' 'OPEN' 'RESV' 10 20 40 /
/
WCONPROD
'PROD1' 'OPEN' 'RESV' 999 3* 123 100 /
'PROD2' 'OPEN' 'RESV' 999 3* 123 100 /
/
TSTEP
3 /
/
END

30
tests/wells_manager_data.data
View File

@ -16,6 +16,10 @@ DYV
DZV
10.0 20.0 30.0 10.0 5.0 /
DEPTHZ
121*2000 /
SCHEDULE
WELSPECS
@ -30,11 +34,11 @@ COMPDAT
WCONPROD
'PROD1' 'OPEN' 'ORAT' 20000 4* 1000 /
/
/
WCONINJE
'INJ1' 'GAS' 'OPEN' 'RATE' 100 200 400 /
/
/
DATES
@ -49,10 +53,26 @@ WCONINJE
'INJ1' 'WATER' 'OPEN' 'RESV' 10 20 40 /
/
END
TSTEP
14.0
/
14.0 /
/
WELSPECS
'TEST1' 'G1' 1 1 8335 'GAS' /
'TEST2' 'G2' 10 10 8400 'OIL' /
/
GRUPTREE
'G1' 'SHIP' /
'G2' 'SHIP' /
/
TSTEP
3 /
/
END

18
tests/wells_manager_data_expanded.data
View File

@ -16,6 +16,12 @@ DYV
DZV
10.0 20.0 30.0 10.0 5.0 /
-- The DEPTHZ keyword is only here to satisfy the old parser; content might
-- completely bogus.
DEPTHZ
121*2000 /
SCHEDULE
WELSPECS
@ -66,7 +72,15 @@ WCONINJE
TSTEP
14.0
/
14.0 /
/
WELOPEN
'INJ1' 'SHUT' 5* /
/
TSTEP
14.0 /
/
END

42
tests/wells_manager_data_wellSTOP.data Executable file
View File

@ -0,0 +1,42 @@
OIL
GAS
WATER
DIMENS
10 10 5 /
GRID
DXV
10*1000.0 /
DYV
10*1000.0 /
DZV
10.0 20.0 30.0 10.0 5.0 /
DEPTHZ
121*2000
/
SCHEDULE
WELSPECS
'INJ1' 'G' 1 1 8335 'GAS' /
'PROD1' 'G' 10 10 8400 'OIL' /
/
COMPDAT
'INJ1' 1 1 1 1 'OPEN' 1 10.6092 0.5 /
'PROD1' 10 3 3 3 'OPEN' 0 10.6092 0.5 /
/
WELOPEN
'INJ1' 'STOP' 5* /
/
TSTEP
10 /
END